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ABSTRACT 


The  three-dimensional  linearized  vorticity  transport 
equations  for  plane  and  pipe  Poiseuille  flow  were  studied 
using  a  highly  generalized  complex  exponential  form  of 
solution  in  both  space  and  time.  The  stability  of  these 
flows  was  examined  using  frames  of  reference  which  move  with 
the  fluid  particles. 
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the  critical  Reynolds  nuabar  is  lowered  by  the  introduction 
of  streamwise  spatial  decay.  This  result  provides  a  new 
basis  for  improving  the  agreement  between  theory  and 
experiment.  Numerical  results  for  pipe  flow  were  not 
obtained  due  to  a  probable  error  in  some  detail  of  the 
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I.   BACKGROUND 

The  problem  of  solving  the  Navier-Stokes  equation  for 
the  onset  cf  instabilities  in  laminar  flow  for  simple 
geometries  has  been  studied  for  many  years.  The  problem  is 
reduced  to  its  simplest  form  by  linearizing  the  vorticity 
transport  eguation,  assuming  sinusoidal  perturbations  in 
space  and  complex  exponential  perturbations  in  time,  and 
reducing  the  eguations  to  their  two-dimensional  form.  The 
solutions  to  this  simplified  problem  are  too  stable  to  agree 
with  experimentally  obtained  values  of  the  critical  Reynolds 
number. 

To  compare  Reynolds  numbers  from  dissimilar  flows,  the 
hydraulic  diameter  and  mean  flow  are  generally  used  as  the 
characteristic  length  and  velocity.  Using  these  criteria, 
the  critical  Reynolds  number  for  a  pipe  is  2600  [Schlichting 
1968]  and  for  a  rectangular  channel  with  a  width  to  height 
ratio  of  8:1  is  also  2600  [ Kao  and  Park  1969].  The  analysis 
in  this  paper  uses  the  channel  semiheight  and  the  pipe 
radius  for  analytic  convenience  (with  the  mean  velocity) . 
Based  on  these  characteristic  lengths  the  experimentally 
obtained  critical  Reynolds  numbers  become  1300  for  a  pipe 
and  730  for  the  8:1  channel.  According  to  Kao  and  Park  the 
critical  Reynolds  number  for  plane  Poiseuille  flow  will  be 
greater  than  that  for  a  channel  of  finite  width. 

The  possible  destabilizing  effect  of  three-dimensional 
flow  in  the  case  of  plane  Poiseuille  flow  was  ruled  out  by 
Sguire's  theorem  [Squire  1933]  which  proved  that 
three-dimensional  perturbations  become  more  unstable  as  they 
are  rotated  into  the  plane  of  the  two-dimensional  flow. 
Generalizing  the  problem  by  allowing  exponential  growth  of 
the  perturbations  in  space  as  well  as  time  greatly  adds  to 
the  complexity,  and  the  solution  of  the  general  problem  has 
not  been  previously  approached. 
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Much  work  has  been  done  for  the  case  of  pipe  flow. 
Flows  with  various  angular  wave  numbers,  and  with  sinusoidal 
spatial  perturhaticns  were  shown  to  be  stable  [Salwen  and 
Grosch  1S72],  perturbations  with  exponential  growth  in 
space,  but  with  strictly  sinusoidal  time  variation  were 
studied  and  found  stable  [Garg  and  Rouleau  1971]  and  the 
combination  cf  exponential  growth  in  space  and  in  time  was 
studied  using  power  series  analysis  [Gill  1973]  and  all  were 
found  to  be  stable. 

Because  of  the  general  agreement  that  pipe  flow  is 
stable  to  inf icitesiraal  disturbances,  two  other  approaches 
to  the  problem  of  accounting  for  the  critical  Reynolds 
number  have  been  used.  First,  it  has  been  assumed  that 
perturbations  which  are  stable  while  they  are  infinitesimal, 
become  unstable  when  they  are  finite.  Studies  have 
concluded  that  finite  disturbances  are  unstable  for  pipe 
flow  [Davy  and  Drazin  1968].  Similar  theoretical 
examination  of  plane  flow  has  reached  the  same  conclusion, 
that  is,  finite  perturbations  are  less  stable  than 
infinitesimal  perturbations  [  Mclntire  and  Lin  1971]. 
Second,  it  has  been  postulated  that  the  origin  of  the 
critical  instabilities  occurs  in  the  entrance  region  of  the 
pipe  (or  channel)  and  both  theoretical  and  experimental 
studies  have  been  done  to  demonstrate  this  [Huang  and  Chen 
1973]  [leite  1956]. 

These  two  explanations  for  the  observed  instability  of 
Poiseuille  flows  have  been  postulated  as  a  result  of  the 
inability  cf  the  linear  theory  to  account  for  the 
experimental  facts.  Unfortunately,  a  totally  general 
solution  to  the  linear  problem  has  never  been  accomplished. 
In  this  paper  such  a  solution  is  presented.  Perturbations 
which  have  a  fully  complex  exponential  form  and  which  are 
fully  three-dimensional  are  introduced  and  the  resulting 
equations   are   solved   using    finite-difference    methods. 
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Results  for  channel  flow  indicate  that  by  varying  the  real 
part  of  the  spatial  wave  number,  the  corresponding  critical 
Reynolds  number  can  be  lowered.  This  is  promising  because 
it  might  help  resolve  discrepancies  between  theory  and 
experiment.  Further  study  is  necessary  to  clarify  these 
results. 


15 


II.   THEORY 
A.   EASIC  EQOATIONS 

The  unperturbed  laminar  flow  of  an  incompressible  fluid 
with  constant  viscosity  is  governed  by  the  vorticity 
transport   and   continuity   eguations.    These   involve   the 

velocity   and   vorticity   vector  functions  V  and  £1   which  are 

well  known  for  the  two  simple   geometries   being   considered 

(see   eguations   (2-3)   through   (2-6) ) .    If   the   flow   is 

slightly   perturbed   and   the   perturbation    velocity    and 

vorticity    are   denoted   by   v   and  o»    respectively,   the 

perturbation  continuity  eguation  and  the  linearized 
vorticity  transport  eguation  are  easily  derived  (see 
Appendix  A) .  These  eguations,  in  their  nondimensional  form, 
may  be  written 

v«v  =  0  (2-1) 

and 


1    v2w   +  o>»vV   -    V«vcu  +    v«v£Z 


-  ^vv    -   3w/3t   =    0  (2-2) 

where  u>  =  vxv.  For  the  two  cases  under  consideration,  plane 
Poiseuille  flow  and  pipe  flow,  the  coordinate  axis 
orientation  is   given   in   figures   2-1   and   2-2   with   the 

velocity  and  vorticity  distributions  V  and  £1. 
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For  Plane  Poiseuille  Flow 


V  =  i3(1-y2)  =  iU(y) 


(2-3) 


&(y)  =  vxV  (y)  =  k3y 


(2-4) 


Figure  2-|     Plane   Poiseuille  Flow 


For    Pipe    Flew 


V(r)    =    e   0  (r)    =   e    2(1-r2) 
x  x 


(2-5) 


X2(r)    =   vxV  (r)    =   e   4r 

G 


(2-6) 


Figure  2-2   Cylindrical  Poiseuille  Flo 


w 
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The  vorticity  transport  equation  in  three  dimensions  is 
actually  three  separate  equations.  These  three  equations 
are  not  independent  (as  shown  in  Appendix  B)  and  represent 
only  two  fully  independent  conditions.  The  continuity 
equation  is  a  scalar  equation  and  the  relation  between 
vorticity  and  velocity  is  a  vector  equation.  Thus  there  are 
six  independent  equations  in  six  unknowns.  The  unknowns  are 
the  three  components  of  the  perturbation  velocity  and  the 
three  components  of  the  perturbation  vorticity. 

These   can   be   reduced   to   three   equations   in   three 

unknowns  by  replacinq  ou  in  equations  (2-1)  and  (2-2)  by  vxv. 
A   further   simplification   can   be   made  by  introducinq  the 

velocity  vector  potential,  W,  where 

v  =  vxW  .  (2-7) 

The  continuity  equation   is   satisfied   identically   because 

v«v  =  v  (vxW)  =  0  (2-8) 

by  a  well  known  vector  identity.   If  v  is   now   replaced   by 

vxW,  the  result  is  two  equations  (the  two  independent 
conditions  cf  the  linearized  vector  vorticity  transport 
equation)  in  three  unknowns.  The  unknowns  are  the  three 
components  cf  the   velocity   vector   potential.    From   this 

result   it   may   be   deduced   that  the   components  of  W  are 

redundant.  This  is  indeed  the  case  and,  as  shown  in 
Appendix  C,  one  cf  them  may  be  set  arbitrarily  to  zero. 

A  useful  way  to  take  advantaqe  of  the  linearity  of 
equation  (2-2)  is  by  seekinq  solutions  which  are  complex. 
Both  the  real  and  iraaqinary  parts  of  any   solution   obtained 
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will,  by  themselves,  be  solutions.  By  choosing  solutions  of 
an  exponential  form,  as  in  equations  (2-9),  the  system  of 
partial  differential  equations  is  reduced  to  ordinary 
differential  equations.  The  desired  form  for  solutions  in 
cartesian  coordinates  is 

W(x,y,z,t)  =  £if(y)  +  jg  (y)  +  kh  (y)  ]e  (2-9a) 


X 
v(x,y,z,t)  =  [iu(yj  +  jv  (y)  +  kw  (y)  ]e  (2-9b) 


X 
w(x,y,z,t)  =  [iC(Y)  +  3V(Y)    +Ck  (y)  Je  (2-9c) 

where 

X  =  ax  +  /5z  +  yt      .  (2-10) 

For  the  cylindrical  case,  the  same  equations,  (2-9)  and 
(2-10),  apply  with  x,y,z  replaced  by  x,r,0.  The  wave 
numbers  a,  /3,  and  y  are,  in  general,  complex.  The  functions 
f ,  g,  and  h  are  defined  by  equation  (2-9a) .    They   are   the 

part  of   the   components   of   W   which  contain  the  y  (or  r) 

dependence  and,  in  general,  are  complex.   Similarly,  u,  v  and 

w   are  part  of  the  components  of  u.   And  £,  rj   and  £  are  part 

of  each  component  of  w,  all  complex  functions  of  y  (or  r)  . 

B.   CAETESIAN  COOEEIHATES 

The  plane  flow  case  for  cartesian  coordinates  is 
considered  first.  Using  equations  (2-9)  to  substitute  into 
equation  (A-13)  ,  as  described  in  Appendix  D,  results   in   an 

equation  in  terms  of  H. 
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-1  vxvxvxvxW  +    vxVxvxvxW  -  vxUxvxW 
Tie 


-     7X7X9^  =  0 

HI 


(D-1) 


Appendix  D  then  develops  the  form   of   the   coefficients   of 
eguation  (D-1)  as  matrices.   The  final  eguation  is 


(Df)  (f 


[.D  I 


(D-55) 


Where  D,  D2  ,  ...  are  the  differential  operators  with  respect 
to  y.   The  matrix  coefficients  are  given  in  eguations  (D-46) 
through  (D-54) .   This  vector  eguation  may   be   expressed   in 
abbreviated  form  by 


r  = 


r  >  =  o 
y 


(B-1) 


as   discussed   in   Appendix   B. 
three  separate  eguations 

r    =  0 
x 


Eguation  (3-1)  is  actually 


(B-2a) 


r    =  0 


(B-2b) 


r     =  0 


(B-2c) 


Examination   of   the   matrix   coefficients   of   eguation 
(D-55)   shows   that   the   highest   order   derivatives  of  the 
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components  f,  g,  and  h,  of  W  in  each  of  the  three   equations 
is  as  fellows 


highest      highest      highest 
eguation     order  of  f   oraer  of  g   order  of  h 

f  4  3  2 

x 

r  3  2  3 

y 

r  234 

z 


TABLE  2-1 


As   demonstrated  in  Appendix  C,  one  of  the  components  of 

W  may  be  set  identically  to  zero,  but,  the  choice  of  h  (y) =0 
leads  to  a  degeneracy  in  the  eguations  for  the  classical 
case  of  fi=Q.  Since  it  is  desirable  to  have  the 
three-dimensional  case  reducible  to  the  plane  flow, 
two-dimensional  case,  it  is  stipulated  that  either  f  or  g  be 
set  to  zero.  Appendix  F  develops  the  form  of  the  boundary 
conditions  for  the  three  cases  of  setting  f,  g,  or  h  to 
zero.   lor  g{y)=0,  the  boundary  condition 

f(±1)  =   h(±1)  (F-6a) 

is  more  complicated  than  for  the  case  of  f(y)=0.  It  is 
therefore  desirable  to  set  f(y)  to  zero.  For  f(y)=0,  the 
boundary  conditions  are 

g(±1)  =  0  (F-4a) 

h(±1)  =  0  (F-Ub) 
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Dh(±1)  =  0  .  (F-Uc) 

From  table  2-1,  derivatives  of  f  (y)  occur  to  the  fourth 
order.  Eguations  (F-4b)  and  (F-4c)  give  four  boundary 
conditions   for  f  (y) .   Table  2-1  also  shows  that  g(y)  occurs 

to  the  third  order  in  eguations  T     and  T  /  while   there   are 

x      z 

only   two   boundary   conditions   in   eguation  (F-Ua).   It  is 


necessary  tc  reduce  the  order  of  g  in  F   and  T   .   This   can 

x       z 

be   done   by   taking   a   linear   combination   of  T     and  V   . 

x      z 

Examination  of  the  matrix  of  eguation  (D-47)  shows  that   the 


coefficient  of  D3g  in  eguation  T     is   a.   In  equation  T     the 

x    Re  z 

coefficient  of  D3g  is  P .   Thus,  the   following   combination 

Re 

of  eguations  T  and  T     will  eliminate  the  terms  in  D3g. 
x      z 

-£T  +  aV     =   0  (2-11) 

x     z 

Appendix   B   shows  that  this  linear  combination  of  eguations 

r   and  T  ,  and  the  condition 
x      z 

r  =  0  (2-12) 

y 

are  sufficient  to  satisfy  the  vorticity   transport   eguation 
written  as  in  eguation  (B-1) ,  under  the  condition  that 

a2  +  £2  *    0  .  (B-9) 


The   resulting   eguations   to   be   solved  are  (2-11)  and 
(2-12)  in  the  unknowns  g  (y)  and  h  (y) .   Performing  the  actual 
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calculations  indicated  by  equation  (2-11)  reveals  an 
interesting  and  extremely  convenient  fact.  Not  only  the 
third  order  derivatives  of  g  disappear  from  equation  (2-11), 
but  g  and  all  its  derivatives  cancel  out  of  the  equation. 
Setting  f  (y)  to  zero,  eguation  (2-11)  only  contains  h  (y)  . 
Thus,  the  eguaticns  are  uncoupled  and  (2-11)  may  be  solved 
for  h  (y)  alone.  With  h  (y)  known,  equation  (2-12)  may  then 
be  solved  for  g  (y) .   Writing  out  equation  (2-11)  yields 

-  a  D«h  +  (aT- a  (a2+P2))D2h  +  (3a2+ a(a2+/32)  T)  h 
Be  $e     H 

+  y(aD2h  +  a(a2+/32)h)  =  0  (2-13) 

where 

T  =  aU  -  1  (a2  +  /32)  =  3a(1-y2)  -  1  (a2+/32)         (D-54) 
le    H  2  He 

The  boundary  conditions  for  equation  (2-13)  are 

h(±1)  =  0  (F-4b) 

Dh(±1)  =  0  .  (F-4c) 

This  homogeneous  differential  equation  is  solvable  as  an 
eigenvalue  problem.  This  technique  is  discussed  in  the 
section  on  numerical  methods.  The  associated  problem  is  to 
solve  eguation  (2-12)  for  g(y).  The  coefficients  of  (2-12) 
may  be  obtained  directly  from  equations  (D-48)  through 
(D-54) .  This  inhcmogeneous  equation  may  be  solved  by 
various  integration  techniques  to  obtain  the  associated 
function  of  g  (y)  for  each  of  the  eigenvalues  obtained  from 
the  solution  of  equation  (2-13).  Solving  equation  (2-12) 
would  be  cf  interest  if  the  functional  shape  of  the 
perturbation  velocities  and  vorticities  were  desired,  but 
the  problem  was  net  dealt  with  in  this  paper. 
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C.   CYLINDRICAL  COORDINATES 

The  eguations  for  pipe  flow,  in  cylindrical  coordinates, 
are  developed  in  Appendix  E.  Equation  (D-1)  is  still  valid 
and  in  abbreviated  form  may  be  expressed  by 


r  =< 


r  >  =  o 

r 


e^ 


(B-11) 


as   discussed   in   Appendix   B.    This   vector   equation   is 
actually  three  separate  eguations. 

(B-12a) 
(B-12b) 
(B-12c) 


r 

x 

0 

r 

r 

= 

0 

re 

= 

0 

The  final  matrix  form  of  equation   (B-11)   is   given   in 
eguaticn   (E-55)   with   the   coefficient  of  equations  (E-40) 
through  (E-48) .   Examination   of   these   coefficients   shows 
that   the   highest  order  of  derivatives  of  the  components  f, 

g,  and  h,  of  H  in  each  of  the  three  eguations  is  as  follows 
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highest  highest  highest 

equation    order  of  f  order  of  g  order  of  h 

T  4  3  2 

x 


r 

r 


r8 


TABLE  2-2 


As  in  the  cartesian  case,  the  order  of  g(r)  in  equations  T 

x 

and   T   is  too  high  for  the  boundary  conditions.   To  reduce 
e 

the  order  of  derivatives  of  g  (r)  in  T       and  T  a   linear 

x       e 

combination   is   taken.   Checking  the  coefficients  of  D3g  in 

both  eguations  it  is  found  that  the  combination 

-/3T  +  aT  =  0  (2-14) 

r  x     9 

eliminates  the  third  order  derivative  of  g (r)  .  This 
combination  of  eguations  does  not,  in  general,  uncouple  the 
problem  by  eliminating  one  of  the  components  of  the  velocity 
vector   potential   from   the  equation.   The  resulting  matrix 

equation,  with  all  three  components  of  W  (one   of   them   may 

later  be  set  to  zero)  can  be  represented  in  the  same  form  as 
equation  (D-55) . 

(D*f)  (D3f)  (D2q) 

l**Wl   +  [H3]iD3g|  +  ^\^rl\l)\lH\ 

*     ([M^/E^])  JDgJ    ♦     <[Mo]+y[No])  jgj    =0  (D-55) 
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but      now    the    matrix    coefficients   are    2X3    matrices,    where   the 

top    row    is    the    combination    of   equations    given    in    (2-14)       and 

the      bottom      row      is      equation      T  .       The   coefficients    are   as 

r 

follows. 


£MJ    = 


0 

r*Fe 

0 


-  a 
He 


(2-15) 


[M3] 


2$ 

r*Re 

a 

Fe 


-    2a 

rTTe 


rRe 


(2-16) 


[M2]    = 


ft-    ft  -£t 

Re    r  JRe  r 

a 
rle" 


-4 


Me 


-t 
Ee 


aT+    3a   -at 
r^Re  Ke 

2$ 
r^Tie 


(2-17) 


[Mt] 


-_#T-    3fl3  +     ft    +  qgft  4a/3        aT+3a/S2-    3a   -    a 

r7      r^Re  f*Ie  r 2 Re        r3Re        r      r-*Re   r  JRe  f"R~ 

r      r^Re 


•aT-     a 

r^Re 


Bl   -   a2 
rJRe   r"Re 


i"0] 


-fttl  +    4g3      2a#T-4a/3  -2a/3t 
r        rbKe        r2      rTRe   r^Re 


4/92-2ag2 
r      r3Re 


tT+  a2    -   £2 
r*Re   r^lfe 


(2-18) 


a(t-1    )  T+    at    +    3a 
F2        r^Re    r*Re 


-  £T-4a/3+ 


■^ 


^ 


_+2a2£ 
►Re  ir*Ae 


(2-19) 


£*23 


CNi3 


r" 

0      a" 

0 

0      0 

- 

F2 

0      a 

r 

a 

0   -$ 
r 

(2-20) 


(2-21) 
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[NQ]  = 


-fit  2aB      at-1 
r     r"2     r2 


0 


-B 


where 


and 


T  = 


t  = 


aU  -  t 
lie 


a2  +  &l 
r? 


(2-22) 


2a(1-r2) 


"  1  (a2+fi|) 

Ke    r2 


(E-48) 


(E-34) 


Since  the  parameter  /3,  which  determines  the  angular 
periodicity,  only  has  discrete  values  as  determined  by  the 
boundary  conditions  in  Appendix  G, 


£  =  ni#    n=0, 1,2, . . . 


(G-4) 


each  value  of  /3  is  considered  independently.  For/3=0,  it  is 
found  from  the  coefficients  (2-15)  through  (2-22)  that  the 
eguations  uncouple,  even  more  completely  than  in  the 
cartesian  coordinate  case.  In  the  first  eguation 
represented  by  these  coefficients  (the  top  row  of  the 
matrices),  both  g  (r)  and  f(r)  drop  out.  The  coefficients 
remaining  represent  a  fourth-order  ordinary  homogeneous 
differential  eguation  in  h (r) .  The  boundary  conditions  for 
this  problem  are  derived  in  Appendices  F  and  G.  Taking  f  (r) 
to  be  identically  zero,  the  boundary  conditions  at  the  wall 
are 


h(1)  =  0 


(F-15b) 


Dh(1)  =  0 


(F-15c) 


The  boundary  condition  at  r  =  0  from  (G-'lO)  is  that  all  even 
derivatives  of  h  are  zero.  In  the  second  eguation  h (r) 
drops  out  because  all  coefficients  contain  /3.  With  f  set 
identically   to   zero,   this  eguation  becomes  a  homogeneous, 
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second  order  differential  equation  in  g (r)  .  The  boundary 
condition  at  r= 1  is 

g(1)  =0  (F-15a) 

The  boundary  condition  at  r=0  from  (G-40)  is  that  g(0)  and 
all  even  derivatives  of  g  are  zero  at  the  origin. 

The  consequence  of  the  complete  separation  of  the 
components  g  and  h  into  two  separate  homogeneous 
differential  equations  is  that  there  are  two  sets  of 
eigenvalues  which  are  derived  independently  of  one  another. 
The  velocity  in  terms  of  f,  g,  and  h  is  derived  in  Appendix 
F.   with  /3=0  and  f(r)=0,  equations  (F-13)  become 

u(r)  =  JD  (rh  (r))  (2-23a) 

r 

v(r)  =  ah(r)  (2-23b) 

w (r)  =  ag  (r)  (2-23c) 

Thus  the  axial  and  radial  components  are  related  to  the 
function  h (r) ,  while  the  angular  velocity  is  completely 
independent . 

Both  of  these  homogeneous  differential  equations  can  be 
solved  as  eigenvalue  problems.  This  technique  is  discussed 
in  the  section  en  numerical  methods.  From  the  eigenvectors 
of  f  and  g  the  velocity  and  vorticity  perturbation  forms  may 
be  found.   This  was  not  done  in  this  paper. 


If  n=1,  that  is,  if  B  =    ft     +  i,  then  the   two   equations 

R 

represented   by   equations   (2-15)  to  (2-22)  do  not  uncouple 

and  must  be  solved  simultaneously,  with  two  of  the 
components  of  the  velocity  vector  potential  as  the  unknowns. 
From  Appendix  F  it  may  be  observed  that  the  mast  complicated 


28 


boundary  condition  at  the  wall 

$f  (1)  =  ah  (1)  (F-17a) 

arises  from  setting  g (r)  equal  to  zero.  From  Appendix  G  it 
is  also  found  that  one  of  the  boundary  conditions  for  n=1  at 
r=0  is 

g(0)  =  ih(0)  (G-41) 

These  inhcmogeneous  boundary  conditions  may  be  avoided  by 
the  choice  of  setting  h(r)=0. 

The  result  is  that  the  first  and  second  columns  of  the 
matrices  of  equations  (2-15)  through  (2-22)  give  the 
coefficients  of  f  and  g  for  each  of  the  two  equations  (the 
top  and  bottom  rows  of  the  matrices)  .  The  boundary 
conditions  at  the  walls  are 

g(1)  =  0  (F-19a) 

f(1)  =  0  (F-19b) 

Df(1)  =  0  (F-19c) 

The  boundary  conditions  at  r=0  are  that  f(0)  and  all  even 
derivatives  cf  f  are  zero,  and  that  g(0)  and  all  odd 
derivatives  cf  g  are  zero. 

Thus,  the  problem  for  n= 1  becomes  two  simultaneous 
homogeneous  differential  eguations  which  must  be  solved 
together.  This  problem  is  discussed  in  the  section  on 
numerical  methods. 

D.   STAEILITY  CRITEEION 

Previous  investigations  into  the  question  of  the 
stability  of  Poiseuille  flows  have  used  variour  criteria  for 
determining   the   stability   of   the  flow  from  the  solutions 
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obtained.   In  particular,  the  real  part  of  the  eigenvalue  y 

R 

is   usually   used   when  there  is  no  real-exponential  spatial 

variation   [Salwen   and   Grosch   1972].    Where   exponential 

growth   in   space   has   been  a  part  of  the  problem  [Garg  and 

Rculeau  1S71]  the  real  part  of  the  spatial  wave   number   has 

been   taken   to   give   the  instability.   In  other  cases,  the 

phase  velocity  has  been  used   to  determine   the   stability. 

Clearly,   a   measure   of  the  stability  of  the  flow  must  take 

into  account  both  the  exponential  rate  of  growth   with   time 

and  with  space. 

Since  it  is  the  stability  of  the  perturbed  fluid 
particles  that  is  of  concern,  it  is  natural  to  consider  how 
the  rate  of  growth  or  decay  of  the  perturbations  is  seen 
from  a  coordinate  frame  moving  with  a  particular  fluid 
particle.  For  cartesian  coordinates,  the  fluid  particles 
can  have  velocities  ranging  from  0  to  1.5  depending  on  their 
distance  from  the  walls,  that  is,  depending  on  the  value  of 
y.  The  distribution  of  fluid  velocities  is  the  Poiseuille 
distribution  of  eguation  (2-3) . 

0(1)  =  3(1-y2)  (2-24) 

These  moving  axes  will  have  a  velocity  with  the  mean  flow  in 
the  x  direction.  Let  x*#  y,  z,  t  be  the  coordinates  and  a, 
$,  y*  the  complex  wave  numbers  with  respect  to  the  moving 
axes.  The  form  of  the  perturbation  vector  potential  for  a 
given  eigenvalue  obtained  as  a  solution  is 

H  =  (D9(y)  +  kh(y))exp(ctx  +  /3z  +  yt)  .  (2-25) 

The  complex  frequency  /*  seen  from  this  moving  reference 
frame  will  be  different  than  /  from  the  fixed  frame.  To 
establish  the  relation  between  /•  and  y,  the  perturbation  is 
written  in  the  moving  frame  and  then  transformed  into  the 
form  which  corresponds  to  the  fixed  frame. 
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W    =     (jg    +    kh)exp(ax'    ♦  $z    ♦  y't) 

=     (J9    ♦    kh)exp[ a(x-Ut)  +  /?z    +    y't] 

=     (jg    +    kh)€Xp[ax    +   (3z  +    (y'-aU)t] 

=     (jg    +    kh)exp[ax    +  /3z  +  yt  ]  (2-26) 

Thus 

y'    -   aU    =     Y  (2-27) 

Solving   for   y'  and  splitting  into  real  and  imaginary  parts 
gives 

y»  =  y  +  a  U  (2-28) 

R     E     B 


y*    =    Y     +  a  U  (2-29) 

'i     I     I 

If   y1   is   positive,  zero,  or  negative,  the  perturbation  is 
R 

said  to   be   unstable,   neutral,   or   stable,   respectively. 

Hence,   the   value   of   y'   is   taken  to  be  a  measure  of  the 

R 

stability  of  each  eigenvalue  obtained. 
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III.  NUMERICAL  METHODS 

A.   CARTESIAN  COORDINATES 

The  eguation  to  be  solved  is  the  linear   combination  of 

the   first   and   third  components  of  the  vorticity  transport 

equation  which  uncouples  the  dependence  of   h (y)   and   g (y)  . 

This   is   developed   in   Section  II  and  results  in  eguation 

(2-13)  . 

-a_D*h  +  (aT-a_(a2+i82))D2h  +  (3a2+  a(a2+£2)  t)  h 
ii€  Pe 

+  y(aD2h  +  a(a2+£2)h)  =  0  (2-13) 

where  T  is  defined  by  eguation  (D-54)  . 

T  =  au  -  1  (a2+£2)  =  3a(1-y2)  -  1  (a2+/32)         (D-54) 
ITe  "2  le 

The  boundary  conditions  are  derived  in  Appendix  F. 

h(±1)  =  0  (F-Ub) 

Dh(±1)  =  0  .  (F-4c) 

Various  methods  are  applicable  to  the  solution  of  this 
homogeneous  fourth-order  eguation.  The  method  of  solution 
applied  is  the  method  of  finite  differences,  approximating 
the  function  h  (y)  by  n  discrete  evenly  spaced  unknowns,  and 
solving  fcr  the  eigenvalues  and  eigenvectors  of  the  system 
of  eguaticns  thus  generated. 

One  of  the  variables,  a,  $  or  y  must  be  chosen  to  be  the 
unknown  eigenvalue.  Since  y  is  linear  while  a  and  0  occur 
to  higher  powers,  the  choice  of  y  produces  a  more  standard 
problem.  Thus,  y  is  the  unknown  and  values  must  be  assigned 
to  a  and  yS.  Equation  (2-13)  has  already  been  written  in  the 
appropriate  form  for  this  method  of  solution,  with  the 
factors  of  y    separated  from  the  rest  of  the   eguation.    In 
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an   abbreviated,   but   general,  form,  equation  (2-13)  can  be 
written 


C  D*h  +  C  D2h  ♦  C  h  +  X(C«D2h  +  C'h)  =  0 

4  2  0  2  0 


(3-1) 


where  the  coefficients  (C's)  are*  in  general,   functions   of 
a,  P,  and  y. 
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Figure  3-1  Finite  Difference  Mesh 

The  range  of  y,  from  +1  to  -1  across  the  channel  height, 
is  divided  into  a  one-dimensional  computational  mesh  of  even 
spacings.  As  shown  in  Figure  3-1,  there  are  n  interior 
points,  n+ 1  divisions  and  n  +  2  total  mesh  points,  including 
the  boundaries.   The  grid  spacing  is  Sy. 

If   the   values  for  the  function  h (y)  at  the  grid  points 

i-2,  i-1,  i+1  and  i+2  are  expanded  as  Taylor   series   about 

th 
the   i     grid   point,   the   resulting   set   of  simultaneous 
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equations  may  be  solved  for  the  second-order  central 
difference  approximations  of  the  derivatives.  These  are  as 
follows. 


Dh.  =  (h.   -h.   )/2(8y)  (3-2a) 

1      x+1   1-1 


D2h   =  (h    -2h  +h    )/(8y)2  (3-2b) 

i      i  + 1    i   i- 1 


D3h.  =  (h.   -2(h.   -h.   )-h.  o)/2(8y)3  (3-2c) 

l      1+2     1+1   1-1    1-2 


D*h   =  (h    -4h    +6h  -4h    +h    )/(Sy)*  (3-2d) 

i     i+2    i+1    i    i-1   i-2 

The   error   cf   eguations  (3-2)  is  of  the  order  of  magnitude 
(By) z. 


The   values  of  h  at  the  boundaries  h(±1)  are  known  to  be 

zero  by  the  boundary  condition  of   equation   (F-4b) .    Thus, 

there   are   n  unkncv:ns--the  values  of  h  (y)  at  the  n  interior 

points — denoted  h  ,  h  ,  ...,   h  ,  ...,  h  .   By   substituting 

12         i        n 

eguations  (3-2)  for  the  derivatives  of  h  in  equation  (2-13) , 

the  fourth  crder  differential  eguation  becomes  a   set   of   n 

linear,   algebraic,  difference  equations  which  form  a  banded 

matrix  with  a  band  width  of  five.    Each   equation,   written 

th 
for  the  i    mesh  point,  includes  the  values  of  h  at  the  mesh 

points  i-2,  i-1,  i,  i+1,  and  i+2.   There   are   four   special 

cases   which   must  be  handled  using  the  boundary  conditions. 

The  eguations  written  for  the  points  2  and  n- 1   include  the 

boundaries.    Since   the   values   for  h  at  the  boundaries  is 

known  to  te  zero,  those  terms  drop  out  of   those   equations. 

The   equations   written  for  points  1  and  n  contain  terms  for 

the  boundaries  and  terms   for   virtual   points   outside   the 

boundaries,   which,   in  accordance  with  the  labeling  scheme. 
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may  be  called  the  -1  and  n+2  values  of  h  (see  Figure  3-2). 
The  second  boundary  condition,  applied  at  the  zeroth 
(boundary)  point,  yields  the  condition 


h    =  h 
-1     1 


(3-3) 


Thus,  in   the   difference   equation   written  at   the   first 

interior  point,  the  term  for  the  boundary  is  dropped  and  the 

coefficient  of   the   virtual   point   is   combined   with   the 

coefficiert   at   the   first   interior   point.    The  equation 

th 
written  for  the  n    point  is  treated  similarly. 

The  matrix  thus  formed  is  broken  into  two  matrices,  one 
containing  the  coefficients  of  the  eigenvalue  y,  and  the 
other  containing  the  remaining  terms.  The  resulting  matrix 
eguaticn  may  be  written 


£x][vj   -  rm{v]  =  o 


(3-4) 


Matrix  [Y]  will  be  a  three-banded  matrix  because  y  only 
occurs  as  a  factor  of  D2h  and  h.  jvj  is  the  eigenvector  of 
values  of  the  function  h  at  the  mesh  points.  Figure  3-3 
shows  the  structure  of  the  matrices. 
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Figure  3-3    Matrix  Form  of  Equation 
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The  two  subprograms,  FUNCTION  CHM1E1  and  FUNCTION  CHM2E1 
supply  the  coefficients  of  these  matrices.  The  input 
parameters,  K  and  Y  are,  respectively,  the  relative  position 
in  the  band  for  each  row,  and  the  value  of  y  for  the 
difference  equation.  The  subprograms  contain  function 
statements  which  pioduce  the  coefficients  obtained  from 
equation  (2-13).  The  parameter  K  determines  how  those 
coefficients  are  combined  to  yield  the  coefficients  of  the 
points  i-2,  i- 1 ,  i,  i+1,  and  i+2,  which  are  determined  using 
the  central  difference  equations,  (3-2) .  Because  the 
standard  eigenvalue  equation  has  a  minus  sign,  as  in 
eguation  (3-4),  the  signs  of  the  coefficients  of  CHM1E1  are 
reversed  from  the  sign  of  the  terms  in  equation  (2-13). 

The  eigenvectors  obtained  from  this  problem  have  either 
even  or  odd  symmetry,  that  is,  they  have  symmetry  about  y=0 
such  that  either 

My)  =  h(-y)  (3-7a) 

or 

h(y)  =  -h(-y)  .  (3-7b) 

Using   this   fact,   the   problem  size  can  be  cut  in  half  by 

solving  for  the   values   of   h    across   only   half   of   the 

i 

channel,   once  with  boundary  conditions  at  y=0  corresponding 

to  even  symmetry  and  a  second  time  with  boundary  conditions 
corresponding  to  odd  symmetry.  This  gives  the  full  set:  of 
eigenvalues  and  vectors.  The  boundary  conditions  for  the 
mesh  points  near  the  center  of  the  channel  are  easily 
derived  fiom  the  conditions  in  eguation  (3-7) . 

Subroutine  MSET  sets  up  the  coefficients  of  either 
matrix  X  or  matrix  Y.  The  variable  MODE  is  used  to  control 
how   the   matrices   are  set  up.   If  MODE=0,  the  matrices  for 
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the  full  channel  result.  M0DE=1  or  MODE=2  will  result  in 
matrices  set  up  for  the  half  channel  with  boundary 
conditions  for  either  odd  or  even  symmetry.  For  the 
half-channel  solution,  the  total  number  of  divisions  in  the 
channel  is  always  even  so  that  the  center  of  the  channel, 
y=0,  falls  between  two  mesh  points.  The  center  of  the 
channel  could  fall  directly  on  a  mesh  point  but  this  would 
result  in  the  inconvenience  of  having  one  more  eigenvalue  in 
the  solution  for  even  vectors  than  in  the  solution  for  odd 
vectors. 

At  the  time  of  writing  there  were  no  computer  routines 
available  to  solve  an  eigenvalue  problem  with  complex 
matrices  in  the  fcrm  of  eguation  (3-4) .  Therefore,  the 
problem  is  ccnverted  to  the  more  conventional  form 

(£Z]  -  y[i])  v  =  o  (3-8) 

by  inverting  matrix  Y  and  forming  the  product 

[Z]  =  £Y][X]~  (3-9) 


The  subroutine  DEIGEO  solves  the  matrix  problem  of 
eguation  (3-4) .  Subroutine  MSET  is  called  to  set  up  matrix 
X.  CDMTIN  is  then  called  to  invert  [X].  CDMTIN  was 
obtained  by  modifying  the  IBM  library  routine  CMTRIN  to  make 
it  applicable  to  double  precision  matrices.  Next,  BMSET  is 
called  to  set  up  matrix  Y  using  space-saving  band  storage. 
The  two  matrices  are  multiplied  using  subroutine  MULDBM. 
Since  all  the  programs  available  for  solving  the  resulting 
eguation  reguire  that  the  real  and  imaginary  part  of  the 
matrix  be  presented  in  separate  matrices,  subroutine  DSPLIT 
is  called  to  accomplish  this.  The  programs  used  to  solve 
for  the  eigenvalues  are  EHESSC  and  ELRH1C  which  are 
available  through  the  International  Mathematical  and 
Statistical  Library.   These  subroutines   reduce   the   matrix 
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into  the  complex  upper  Hessenberg  form  and  then  solve  for 
the  eigenvalues.  Another  program  which  can  be  used  to  solve 
for  the  eigenvalues  is  EISPAC  (Eigensystem  Subroutine 
Package)  which  is  available  from  the  IBM  library.  This 
program  was  also  used  to  solve  this  problem  and  gave  results 
identical  to  the  results  obtained  using  the  I.K.S.L. 
routines.  DEIGEO  solves  the  problem  twice,  once  for  the 
even  eigenvectors  and  once  for  the  odd  eigenvectors.  The 
results  are  then  passed  back  to  the  calling  program. 

Three  main  programs  have  been  written.  Program  #1 
takes  values  of  a  and  )9  as  input,  calls  DEIGEO,  and  outputs 
a  list  of  the  eigenvalues,  along  with  the  stability  for  each 
eigenvalue,  determined  using  equation  (2-35).  A  printer 
plot  is  also  cutput.  Program  #1A  is  the  same  as  program  #1 
except  that  a  plot  on  the  Calcorap  plotter  is  output  using 
subroutine  DRAW.  This  output  is  used  in  this  paper  to 
present  the  graphs  of  the  eigenvalues.  Program  #2  was  used 
to  determine  the  shape  of  the  stability  curves.  This 
program  performs  five  functions  depending  on  the  value  of 
the  control  parameter  MODE.  The  stability  curves  determined 
are   ail   two-dimensional  plots  of  (X      (the  imaginary  part  of 

a)    versus  the  Reynolds   number,   holding   /5  ,  $    ,      a  ,   the 

I     R     R 

velocity,  and  the  stability  constant. 


For  M0DE=1,  the  stability   for   a   given   set   of   input 
parameters,  a,    yS,   Reynolds  number  (REY)  and  velocity  (VEL) 
is  determined  and   printed.    Equation   (2-35)   is   used  to 
determine  the  stability. 

For  M0DE=2  a  point  on  the  line  of  constant  stability 
with  an  input  value  of  STAB  is  determined  for  a  given  input 
value  of  REY.   Referring  to  Figure  3-Ua,  the   initial   guess 
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has   a   value   of  a   and  REY,  and  corresponds  to  point  1.   A 

second  point  (2)  is  determined  by  adding  DAI  to  the   initial 

value   of   a  .    The   stability   at   both  of  these  points  is 
I 

determined  by  calling  the  subroutine  EIGLS  which   determines 

the  value  of  the  stability  of  the  least  stable  eigenvalue 
for  a  given  set  of  input  parameters,  a,  /?,  REY,  and  VEL.  A 
linear   interpolation   on   a    is   made  to  determine  a  third 

guess  (3)   for   which  the   stability   is   also   calculated. 

Subroutine  EFIT2  is  then  called  to  make  a  second-order 
approximation  of  the  value  of  a   which   corresponds   to   the 

desired   value  of  the  stability  (STAB) .   This  final  value  is 

printed. 


Por   1500.5=3,   the   same   operations  are  performed  as  for 
M0DE=2,  except  that  REY  is  varied   instead   of  a    .    Figure 

3-Ub  shows  this.   The  initial  guess  is  point  1,  DREY  is  used 

to  find  point  2,  a  linear  interpolation  to  find  3  and  a 
second-order  fit  to  determine  a  final  estimate,  which  is 
printed. 


For  M0DE=4,  the  point  of  minimum  REY  on  the  neutral 
stability  curve  is  determined.  Two  points  (1b  and  1c)  are 
determined  from  the  initial  guess  (1a) ,  by  adding  and 
subracting  DAI  from  a   .   Their  stabilities  are  determined  by 

calling  EIGLS  and  a  second  order  fit  determines  the  least 
stable  point  (1)  for  the  initial  value  of  REY.  A  new  value 
of  REY  is  picked  by  adding  or  subracting  DREY  as 
appropriate .  This  results  in  point  2a.  The  points  2b  and 
2c  are  found  by  adding  ±DAI.  The  stability  of  these  three 
points  is  determined,  resulting  in  the   least   stable   point 
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(2)  for  this  value  of  REY.  A  third  point  (3a)  is  determined 
using  a  linear  interpolation  from  points  1  and  2,  two  points 

(3b  and  3c)  are  chosen  by  adding  tDAI  and  the  stabilities  of 
these  points  results  in  the  least  stable  point  (3)  . 
Finally,  a  seccnd-order  fit  is  made  using  points  1 ,  2,  and, 
3  to  find  the  point  corresponding  to  zero  stability.  This 
point  will  be  the  point  of  minimum  Reynolds  number  on  the 
neutral  stability  curve. 

For  M0DI=5,  the  point  of  minimum  stability  is  determined 

for  a  set   of   input   values   a  ,   /3  ,  0   ,   and   VEL.    The 

R    I     R 

procedure   is   similar  to  that  for  M0DE=4.   An  initial  guess 

(1a)  and  two  points  (1b  and  1c)  are  used  to  find  point  1. 
This  procedure  is  repeated  twice  at  values  of  REY  determined 
by  adding  or  subtracting  DREY  from  the  initial  guess. 
Points  1,2,  and  3  are  then  used  to  approximate  the  point  of 
minimum  stability  using  subroutine  DFIT2. 

B.   CYLINDRICAL  COORDINATES 

Programs  to  solve  the  two  cases,  n=0  and  n=1,  were 
written  separately.  The  method  of  finite  differences  is 
used  in  each  case.  All  subroutines  used  in  the  cartesian 
case  are  used  in  the  cylindrical  case  with  the  exceptions  of 
the  functions  returning  the  coefficients  and  the  subroutine 
setting  up  the  matrices. 

Functions  CHM1E1,  CHM2E1,  CHM1E2,  and  CHM2E2  return  the 
values  of  the  coefficients  of  the  velocity  vector  potential 
component  h  for  the  five  adjacent  mesh  points  in  the  central 
difference  representation  of  the  equations.  This  is  done  in 
a  manner  exactly  the  same  as  the  functions  CHM1E1  and  CHM2E1 
in  the  cartesian  case.  The  "M2"  means  those  terms  which 
also  contain  the  eigenvalue  y ,  and  are  used  to  set  up  the 
matrix  Y.   The  "HI"  means  those  remaining  terms  which  do  not 
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contain  y,  and  are  used  to  set  up  the  matrix  X.   "E1"   means 

that  the  terms  come  from  the  linear  combination  of  T     and  F 

x      e 

given  in  equation  (2-14)  and  which  are  written   out   as   the 

top   row   of   the   matrices   of  (2-15)  through  (2-22)  .   "E2" 

means  that  the  terms  come  from  the  equation  T  ,         and   which 

r 

are  the  bottcm  row  of  the  matrices  of  (2-15)  through  (2-22) . 

Similarly,  functions  CGM1E1,  CGM2E1,  CGM1E2,  and  CGM2E2 
return  all  the  values  of  the  finite  difference  coefficients 
of  the  component  g  (r) ,  and  CFM1E1,  CFM2E1,  CFM1E2,  and 
CFM2E2  return  the  values  of  the  coefficients  of  f (r) . 

Subroutine  I1SET2  performs  exactly  the  same  function  as 
MSET,  with  the  exception  that  MSET2  allows  the  matrix  to  be 
positioned  with  a  starting  value  ether  than  (1,1) .  This  is 
necessary  because  for  cases  other  than  n=0,  two  equations 
must  be  sclved  simultaneously  and  matrices  must  be  placed 
adjacent  to  each  other  to  accomplish  this.  Also,  MSET2 
forms  a  mesh  which  includes  the  point  r=0,  in  contrast  to 
MSET  which  had  the  center  of  the  channel  fall  between  mesh 
points.  The  boundary  conditions  must  be  written 
appropriately.  MSET2  sets  the  correct  boundary  conditions 
at  the  wall  r=1,  but  does  not  set  any  boundary  conditions  at 
r=0  since  they  are  different  for  each  case.  Thus,  the 
boundary  conditions  at  r=0  must  be  set  by  the  calling 
program. 

!•   Axisymraetric  Flow  (n-0) 

Prcgram  #R1A  solves  the  first  equation  for  the 
function  h (r) .  As  shown  in  part  II,  section  C,  the  boundary 
conditions  can  be  written 

h(1)  =  0  (3-10) 

Dh(1)  =  0  (3-11) 
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h(0)  =0  (3-12) 

D2h(0)  =  0  (3-13) 

All  even  derivatives  at  the  origin  are  zero,  but  only  (3-12) 
and  (3-13)  are  needed  to  satisfy  the  difference  equations  at 
r=0.  The  above  equations  imply  that  the  boundary  mesh 
points  at  r=0  and  r=1  are  left  out  of  the  difference 
equations.  When  the  virtual  point  outside  the  boundary  r=1 
appears,  its  coefficient  is  added  to  the  first  mesh  point 
inside  the  boundary.  When  the  virtual  point  outside  the 
boundary  r=0  appears,  the  negative  of  its  coefficient  is 
added  to  the  first  mesh  point  inside  the  boundary.  The 
boundary  conditions  at  r=1  are  set  by  subroutine  MSET2  while 
the  boundary  conditions  at  r=0  are  set  by  the  main  proqram 
#R1A.  The  rest  of  the  solution  proceeds  exactly  the  same  as 
the  method  used  for  the  cartesian  case  discussed  in  section 
A  of  part  III. 

Frogram  #B1B  solves  the  second  equation  for  the 
function  g (r) .  The  boundary  conditions,  from  part  II 
section  C,  are 

g(1)  =  0  (3-14) 

g(0)  =0  (3-15) 

In  addition  to  (3-15),  all  odd  derivatives  of  g  are  zero  at 
r=0,  but  cnly  (3-15)  is  required  to  satisfy  the  difference 
equations  at  r=0  since  the  highest  derivative  of  g  is  the 
second  derivative  and  the  result  is  a  tri-diagonal  matrix. 
Again,  the  solution  procedure  is  identical  to  the  cartesian 
case. 

Both  programs  #R1A  and  #R1B  output  a  list  of  the 
eigenvalues  and  the  stability  for  each  eigenvalue.  Printer 
plots  of  the  eigenvalues  are  also  output. 
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Figure   3-5    Finite    Difference   Matrix    for   Pipe    Flow 
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2.   njj 

For  this  case  the  two  equations  must  be  solved 
simultaneously.  This  is  done  by  writing  two  sets  of  n 
equations  for  the  set  of  2n  variables  including  n  unknown 
values  of  f  and  n  unknown  values  of  g.  The  function  h  is 
set  identically  to  zero  as  described  in  section  II,  part  C. 
The  boundary  conditions  are 

f(1)  =  0  (3-16) 

Df(1)  =  0  (3-17) 

g(1)  =  0  (3-18) 

f(0)  =0  (3-19) 

D2f  (0)  =  0  (3-20) 

g(0)  =  0  (3-21) 

The  matrices  are  set  up  as  shown  in  figure  3-5  using 
subroutine  MSEI2  four  times.  The  upper  left  section  is  set 
up  using  function  CFM1E1.  CGM1E1  is  used  for  the  upper 
right  section,  CFM1E2  for  the  lower  left  and  CGM1E2  for  the 
lower  right.  A  second  matrix  is  set  up  using  an  identical 
procedure,  but  for  the  "H2"  coefficients,  that  is,  the 
coefficients  which  contain  /.  The  procedure  is  then  to 
invert  this  seccnd  matrix,  multiply  it  with  the  first  and 
solve  for  the  eigenvalues.  Program  #R2  carries  out  the 
solution  and  then  lists  the  eigenvalues  found,  with  their 
stability.   A  plot  is  also  output. 
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IV .   RESULTS 

A.   COMPARISON  WITH  PREVIOUS  RESULTS 

1 .  Neutral  Stability  Curve 

The  neutral  stability  curve  was   calulated   for   the 

two-dimensional,   sinusoidal  case,  that  is,  for  a   =  /9  = ft  =0 . 

R   I   R 

Figure  4-1  compares  the  results  with   those   of   Grosch   and 

Salwen  (1S68).  Using  n=60  (60  mesh  points  across  the  half 
channel)  gives  almost  the  same  accuracy  as  the  published 
results.  The  computation  time  for  finding  the  eigenvalues 
using  n=60  was  over  one  minute  (on  an  I.B.M.  360),  compared 
to  a  computation  time  of  approximately  10  seconds  for  n=30. 
For  this  reason  a  mesh  size  of  30  was  used  to  obtain  all 
results  presented  in  this  paper  except  where  noted. 

2.  Squire1 s  Theorem 

Sguire*s  theorem  can  be  verified  numerically.  As 
expected,  the  rotation  of  three-dimensional  sinusoidal 
disturbances  onto  the  two-dimensional  plane  always  has  a 
stabilizing  effect.  From  Sguire* s  theorem,  it  follows  that 
the  addition  cf  a  transverse  sinusoidal  component  to  a 
two-dimensional  sinusoidal  wave  is  absolutely  stabilizing  in 
that   the   minimum   Reynolds  number  on  the  neutral  stability 

curve  must  increase  for  increasing  /3  .   Figure   4-12,   which 

I 

plots   the   minimum  Reynolds  number  for  various  values  of  a 

R 

and  B   ,   shows   this   to   be   so.    For   a  =0,   the   neutral 
I  R 

stability  curve  does  not  seem  to  have  another  lobe  (labelled 

e-f-g  in  figure  4-10a),  as  the  curves  do  for  a  *0.   For   the 

R 

cases   for   which  Squire's  theorem  applies,  a  =0,  figure  4-I2 

R 
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shows  that  increasing  /3  is  stabilizing.    Squire's   theorem 

is  not  applicable  to  the  other  cases,  where  a  #0. 

R 


B.   STABILITY  CURVES 

The   stability   criterion   used   is  discussed  in  section 

II-D.   It  is  based  on  how  an  observer  in  a  moving  coordinate 

system   would   experience   the   combination   of   exponential 

growth  in  space  and   in   time   of   the   perturbations.    The 

stability  is  given  by  y1  in  equation  (2-28) . 

fi 

v«  =  v  +  a    u  (2-28) 

B     R     R 

The   stability  observed  from  a  coordinate  system  moving  with 

velocity  U   is   said   to   be  stable,   neutral   or   unstable 

according   to   whether  y*    is  less  than,  equal  to,  or  greater 

R 

than  zero.   The  velocities  of   interest   are  those   of   the 
fluid  particles,  which  vary  from  0  to  1.5. 


The  neutral  stability  curves  for  different  values  of  a  , 

R 

with    ft  =  0,    are    shown    in   figure    4-2.       The   stability    criterion 

I 

used  is  for  a  fixed  coordinate  system  where  U=0  in   equation 

(2-28) .   This  is  the  same  as  the  instability  observed  from  a 

ncnmoving  frame  of  reference.   The  structure  of   the   curves 

for   a  *0   appears   to  be   qualitatively  different  from  the 
R 

curve  for  a   =0  in  that  (referring  to  figure  4-10a)  the   lobe 
B 

e-f-g  is  net  present  on  the  neutral  curves  for  a  =0. 

R 


Figure  4-3  shows  the  effect  on  the  neutral  curves  of  the 


U8 


addition   of  a  /?  component.   It  appears  that  only  the  upper 

lobe  (labelled  b-c-d  in   figure   4- 10a)   of   the   curves   is 

affected.    The   stability  is  y%    from  eguation  (2-28) ,  for  a 

R 

fixed  coordinate  system  (U=0) . 


Cases   involving  nonzero  values  of  /3  were  not  computed. 

R 

This  was  primarily  because  of   a   lack   of   time,   but   also 

because   a   hypothetical   case   of   this   kind  might  well  be 

difficult  if  not  impossible  to  realize  experimentally. 


Curves   of   constant   stability  for  a  =0  and  a.    =-.04  are 

R        R 

shown  in  figures  4-4  and  4-5.   It  is  important  to  note   that 

eguation   (2-28)   can   be  used  (for  nonzero  values  of  a  )  to 

R 

transform  values  of  calculated   stability   as   seen   by   one 

reference   frame,  into  different  values  of  stability  as  seen 

by  another   reference   frame   moving   with  a   different   U. 

Conversely,  if  a  curve  of  constant  stability  is  given  for  one 

reference  frame,  it  is  possible  to   find   another   reference 

frame  which  sees   that  curve  as  having  any  desired  stability, 

such  as  neutral  stability,  for   example.    The   velocity   of 

this   frame   cf   reference   is   meaningful   only  if  it  falls 

within  the  limits  of  the  fluid   particle   velocities,   0   to 

1.5.   Thus,  for  a  =-.04,   the  neutral  curve  for  a  coordinate 

R 

frame  moving  with  the  mean  velocity,  U=1.0,   corresponds   to 

the   curve   of   constant  stability  for  /  =.04  as  seen  from  a 

R 

fixed  coordinate  frame.   Correspondingly,  the  neutral   curve 

for  the  maximum  velocity,  0=1.5,  corresponds  to  the  curve  of 

constant  stability  for  /  =.06  for  U=0. 

R 
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C.   EIGENVALUES 

The  method  of  finite  differences  was  used  to  solve  the 
differential  eguations  (described  in  section  III) .  A  mesh 
of  n  finite  difference  grid  points  was  placed  across  the 
half  chanr.el. 

The  effect  cf  mesh  size  on  the  eigenvalues  is 
demonstrated  in  figures  4-6a,b,c,  and  d.  From  these  graphs 
it  can  be  seen  that  the  effect  of  increasing  the  accuracy  of 
the  soluticn  by  increasing  the  number  of  mesh  points  is  to 
stretch  out  the  eigenvalues.  For  high  enough  n,  the 
eigenvalues  form  a  "Y"  lying  on  its  side  with  the  branches 
near  the  y     axis.   The  other  two  branches,  which   move   away 

from   the  y       axis   with   increasing   n,  seem  to  be  "false" 

eigenvalues  which  can  be  made  ever   increasingly   stable  by 

improving  the  accuracy  of  the  solution.  A  true,  qualitative 
picture  of  the  eigenvalues  is  as  shown  in  figure  4-  10b,  with 
the  line  of  eigenvalues  extending  indefinitely  to  the  left. 
Note  also  that  in  all  the  eigenvalue  plots  shown,  the  value 
of  y  corresponding   to  the  base  of  the  "Y"  formed  by  the 

eigenvalues  is  egual  to  the  value  of  a  . 

The  effect  cf  Reynolds  number  is  shown  by  figures 
4-7a#b#c,d.  It  is  apparent  that  changing  Reynolds  number 
has  an  effect  similar  to  that  of  changing  n.  This  seems  to 
indicate  that,  for  a  given  mesh  size,  lower  Reynolds  numbers 
yield  more  accurate  eigenvalues  than  higher  ones.  This  is 
shown  to  be  true  later  in  this  section. 


The   effect   of   a    on   the  eigenvalues  is  shown  by  the 

R 
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series  of  plots  in  figures  4-7b,  4-8a,  and  4-8b.   The  effect 

of  decreasing  a   is  to  make  the  flow  less  stable.   Note  that 

R 

both  branches  of  the  eigenvalue  plots  become   unstable.    No 

plots  of  positive  a   are  shown,  but  the  results  are  to  shift 

R 

the    eigenvalues    to    decreasing    y  ,    that    is,    stabilizing. 

R 


The   effect   of   a    is  shown  in  figures  4-9a,  4-7b,  and 

I 

4-9b.   Increasing  a   tends  to   increase   the   f reguency { y  ) . 

I  I 

that  is,  shorter  wavelengths  are  associated  with  higher 
f reguencies . 

The  significance  of  the  different  parts  of  the  stability 

curves  is  shown  in  the  series  of  figures,  4-11.   The  neutral 

stability   curve   for   a  =-.04  is  shown  in  figure  4-10a  with 

R 

the  parts  of  the   curve   labelled.    Figure   4- 10b   shows   a 

typical  eigenvalue  plot  for  large  n  with  the  branches 
labelled.  Eigenvalues  for  the  points  labelled  A,  3,  C,  D 
and  E,  are  plotted  in  figures  4-1 1a,b,c,d, e.  From  these 
plots  it  is  observed  that  the  traditional  neutral  curve, 
represented  by  the  instability  of  point  B,  is  associated 
with  branch  LF  of  the  eigenvalues  corresponding  to  lower 
f reguencies.  The  curve  a-b  d-e  that  cuts  vertically  through 
the  upper  lcbe  is  associated  with  the  instabilities  of 
pcints  A  and  D.  These  points  have  instabilities  on  the 
other  branch  (KF)  of  the  eigenvalues,  the  one  associated 
with  higher  frequencies.  Figure  4-1 1c  illustrates  that  at 
point  C,  both  branches  of  the  eigenvalues  are  unstable.  At 
point  E,  figure  4-11e  shows  that  the  two  branches,  HF  and 
LF ,  of  the  eigenvalues  have  almost  disappeared,  which  seems 
to  relate  to  the  sudden  decrease  in  stability  of  the  lower 
lobe  e-f. 
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D.   STABILITY 

The  instability  of  the   top   lobes   (labelled   b-c-d   in 

figure   4-10a)   of  the  neutral  curves  of  figures  4-2  and  4-3 

is  qualitatively   different   from   instabilities   elsewhere. 

For   this   reason,   and   because   these   lobes   are  directly 

analogous  to  the  neutral  curve  for  a  =0,  this  region  of   the 

R 

curves   is   considered   separately  from   the   rest   of   the 

instabilities.   Each  of  these  lobes  has  a  point   of   minimum 

Reynolds   number.    These   are   plotted   in   figure  4-12  for 

different  values  cf  a   and  y9  .   Once   again,   the   stability 

R      I 

criterion  is  based  on  eguation  (2-28)  for  a  fixed  coordinate 
system. 

Each  of  the  upper  lobes  (b-rc-d)  also  has  a  maximum 
instability  associated  with  it  similar  to  the  maximum 
instability  shown  in  figures  4-4  and  4-5.  These  values  are 
shown  in  figure  4-13  which  plots  the  coordinates  (a   vs.  Re) 

of   the   Eaximum   instabilities.    Only   the  case  of  /5  =0  is 

considered.   Note  that  the  values  of  a   plotted  only   go  to 

R 

-.04.   This  is  because  the  instabilities  associated  with  the 


part  of  the  neutral  curve  labelled  a-b  d-e  in   figure   4-10a 
overwhelm   the   instab 

increasingly  negative. 


overwhelm   the   instabilities   of   the  lobe  b-c-d  as  a      gets 

R 


If   the   values  of  the  maximum  instabilities  are  plotted 

vs.  or  ,   as  in  figure  4-14,  they   form   the   curve   labelled 
R 
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0=0.   The  values  cf  stability  of  figure  4-13  correspond  to  a 

ncnmoving  coordinate  system.    Using   eguation   (2-28),   the 

maximum   instabilities   for   nonzero  velocities  are  plotted. 

The  lines  appear  to  be  almost  straight.    Note   that   for   a 

velocity    cf    approximately    .41,   the   line   of   maximum 

instability  is   approximately   horizontal.    For   velocities 

greater   than   .41   the   line  intersects  the  a      axis  at  some 

R 

negative  value  cf  a  .   For  a   below  this  value,  the  flow   is 

B         R 

stable.    Assuming   that  these  lines  are  straight,  values  of 

velocity  less   than   .41  will   make   the   flow   unstable   at 

Reynolds   numbers   as  low  as  desired  by  choice  of  a  suitably 

negative   value   of    a  .     Although    not    yet    verified 

R 

numerically,   it   seems  possible  that  for  velocities  greater 

than  .41  there  may  be  a  critical  value  of  the  Reynolds  number 

below   which  no  instabilities  are  found.   Curves  of  constant 

stability  for  positive  a   must  be  plotted  to  see  if  this   is 

R 

so.         Cases   of    /3  #0    were   not   considered   due   to   lack    of    time. 

I 


The  most  critical  portion  of  the  neutral  curve  occurs  at 

values    of    a     approaching    zero,   that   is,   for   long 

I 

wavelengths.   The  stability  of  points  along  a  line  close   to 

the  Re  axis  (at  a  value  of  a  =.01)  was  determined  for  values 

I 

of  a   ranging  ficu  -.02  to  -10.0.   Eguation  (2-28)  was   used 
R 

to  translate   these  values   of   the   stability   found  for  a 

ncnmoving  coordinate  frame  of  reference  to  the  velocity  of  a 
coordinate  axis  for  which  the  stability  was  exactly  neutral. 
The  resulting  curves  are  plotted  in  figure  4-15.  The  region 
of  instability  lies  below  and  to  the  right  of  the  curves. 
It  can  be  Eieen    frcm  the   graph   that   there       no   minimum 
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Reynolds   number   associated   with   these  instabilities.   In 

particular,  a  value   of   a  =-10.0   seems   to   be   completely 

R 

unstable   for   all    velocities   and    Reynolds   numbers. 


E.  PIPE  ELCW 

The  results  for  the  numerical  solution  of  the  pipe  flow 
eguations  seem  to  indicate  that  there  is  an  error  somewhere 
in  either  the  analysis  or  the  numerical  solution  technique. 
All  three  programs,  R1A,  R1B  and  R2,  produce  eigenvalues 
which  become  increasingly  unstable  for  decreasing  Reynolds 
number.   Further  study  is  necessary  to  resolve  this. 

F.  NUMERICAL  ACCURACY 

Figure  4-16  shows  the  neutral  stability  curves   for   two 

values   of  a        for   a  mesh  with  30  points.   Portions  of  the 
R 

curve  with  n=60  are  also  plotted.   These  graphs  support   the 

observation   made   earlier   that   decreasing  Reynolds  number 
seems  to  be  associated  with  increasingly  accurate  solutions. 

Calculations  were  done  on  an  I.B.M.  360  computer  in  both 
single  and  double  precision,  that  is,  with  computer  word 
lengths  cf  32  and  64  bits  respectively.  Some  calculations 
were  also  dene  on  an  X.D.S.  9300  with  a  48  bit  word  length. 
With  n  =  60,  the  eigenvalues  obtained  on  the  9300  and  in 
double  precision  en  the  360  agreed  to  four  to  six  places. 
Those  obtained  in  single  precision  on  the  360  sometimes 
agreed  to  only  the  first  place  with  the  double  precision 
values.  Thus  it  appears  that  the  64  or  48  bit  word  length 
is  sufficient  for  matrices  of  size  n=60  while  in  single 
precision  en  the  I.B.M.  360,  some  eioonvalues  have 
significant  errors. 
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Figure  4-14    Values  of  maximum  instability  of  the 

upper   lobe    vs.    Qf_    for   various  velocities 

R 


71 


_  o 


ro 

i 

O 

mmm 

X 

*■"* 

i_ 

0> 

o 

q: 

v> 

JO 

o 

cc 

o> 

a 

£ 

o 

TO 

— 

c 

0) 

D 

x: 

</> 

*~ 

<U 

•*- 

'♦_ 

o 

o 

>^ 

o 

*— 

<u 

■^ 

> 

X> 

o 

to 

3 

(A 

o 

— 

v. 

o 

a 

V- 

> 

3 

a> 

c 

0) 

a: 


i£) 

i 

V. 

3 

en 


(D)      ^IP0|8A 


72 


o 

* 

O    DC 
CO 


CO 


CD 


OJ 


> 

O 


o 

o 

5 


o 

D 

i_ 

U 
O 

o 

o 

u 

E 

3 


CD 

i 

0) 

l_ 

3 


I L 


CM 


_J L 

00 


CD 


73 


V.    CONCLUSIONS  AND  RECOMMENDATIONS 

The  solution  of  the  linearized  vorticity  transport 
eguation  was  simplified  considerably  by  the  introduction  of 
the  velocity  vector  potential.  This  method  has  not  been 
mentioned  in  previous  publications,  at  least  not  in  those 
with  which  the  author  is  familiar. 

The  effect  of  the  number  of  finite  difference  mesh 
points  on  the  eigenvalues  is  important  for  understanding  the 
effects  of  the  ether  parameters.  This  is  explained  in  the 
results  section,  and  is  believed  to  be  information  not 
previously  published. 

The   results   presented   here   show   that   the   critical 

.Reynolds  number  for  plane  Poiseuille  flow  can  be  lowered   as 

much  as  desired  by  proper  selection  of  ether  parameters.   In 

particular,  it  is  highly  significant  that  a   negative   value 

of   a  ,   that  is,  streamwise  decay  in  space,  produces  highly 
R 

unstable  rates  cf  growth   in   time.    The   lowering   of   the 

critical   Reynolds   number   in  this  way  provides  a  new  basis 
for  improving  the  agreement  between  theory  and  experiment. 

Two    separate    modes    of    instability    were    found 

corresponding  to  the  two  branches  of  the   eigenvalue   plots. 

The  classical  case  of  instability  found  for  a  =0  corresponds 

R 

to  only  one  cf   these   branches   becoming   unstable.   It   is 

significant   that   for  each  different  value  of  a  ,  this  mode 

R 

has  a  maximum  instability  associated  with  it.   For   a   fluid 
particle    velocity    of    approximately   .4  1,  this   maximum 


7U 


instability  seems  to  be  invariant  with  a  .    An   additional 

B 

mode   of   instability,   corresponding  to  longer  wavelengths, 

appears  for  nonzero  a    .   No  maximum  value  of  instability   is 

R 

present   as   in  the  first  case.   This  mode  has  been  shown  to 

be   unstable   at   all   values   of   Reynolds   number   for 

sufficiently  negative  values  of  a  . 

R 


The  stability  of  a  fluid   particle   has   been   shown   to 

depend   upon  its  velocity.   The  fact  that  negative  a   yields 

R 

the  greater  instabilities   indicates   that   fluid   particles 

with  the  lowest  velocities,  that  is,  those  nearest  the 
walls,  will  be  the  most  unstable.  This  seems  to  agree  with 
experiment.  For  example,  from  Schlichting  (1968),  the 
transition  is  "characterized  by  an  amplification  of  the 
initial  disturbances  and  by  the  appearance  of 
self-sustaining  turbulent  flashes  which  emanate  from  fluid 
layers  near  the  wall  along  the  tube".  Further  comparison 
with  experiment  is  recommended. 


The  eigenvectors  of  the  solutions  to  the  vorticity 
transport  eguaticn  were  not  investigated.  Those 
corresponding  to  the  least  stable  eigenvalues  should  be 
studied. 

Since  satisfactory  results  were  not  obtained  for  pipe 
flow,  apparently  because  of  an  error,  an  independent  study 
of  this  case  should  be  made,  including  the  boundary 
conditions  and  the  numerical  methods.  It  seems  likely  that 
the  solutions  will  be  qualitatively  similar  to  those 
obtained   for   plane   flow.     If    this    be    true,    then 


75 


instabilities   should  be   found   for  negative  values  of  a  . 

E 

This  would  be  significant  because   it   is   generally   agreed 

that    for   a  =0,   pipe   flow   is  stable   to   infinitesimal 
R 

perturbations. 
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APPENDIX  A 

DERIVATION  OF 
VORTICITY  TRANSPORT  EQUATION 

The   flow   of   an  incompressible   fluid   with   constant 
viscosity  is  governed  by  the  Navier-Stokes  equation 

w2v  -  (v  «v)v  -  vp  /p     -  av  /at  =  o         (A-1) 

d     d    d     d        d   d 

where   the   subscript  d  indicates  the  dimensional  quantities 

of  the  velocity  vector,  pressure,  and  time  (V  , p  ,  and   t  ). 

d  d        d 

The  kinematic  viscosity  v   and  the  density  p    are  constants. 

This  equation  is  nondiraensionalized   using  a   reference 

length   L   equal   to   the  half  width  of  the  channel  in  plane 

flow  or  the  radius  of  the  pipe  in  pipe  flow.   The   reference 

velocity   is   the  averaqe  velocity  of  the  laminar  flow  V 

avg 

The  resulting  eguation  is 


1_V2V  -  (V«v)V  -  vp  -  av/3t  =  0  (A-2) 

T3e 

where  Re  =    IV   / y   is  the  Reynolds   number.    This   flow   is 
avg 

also  governed  by  the  continuity  eguation. 


V«V  =  0  (A-3) 

Other   than   equation   (A-1),  all  the  equations  presented  in 
this  paper  are  in  ncndimensional  form. 

The  vorticity  transport  equation  is  obtained  by  taking 
the  curl  (vx)  of  the  Navier-Stokes  equation.  The  definition 
of  the  vorticity  is 
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ft  f  vxV  (A-4) 

and  the  following  vector  identities  are  applicable; 

vxvp  =    0  (A-5a) 


vX[(V«v)V]  =  vx£v(V«V)/2-Vx  (vxV)  ]  =  -v*(Vxft) 


(A-5b) 


vx(v2V)  =  v2ft 


(A-5c) 


vx(3v/<3t)    =  aft/ at 


(A-5d) 


Combining  equations  (A-4) ,     (A-5) ,  and  the  curl  of  (A- 2) ,  the 
result  is 


1  *2ft  +  7X  (vxft)  -  aft/3t  =  o  . 

ITe 


(A-6) 


Next,  it  is  assumed  that  the  flow  is  made  up  of  the 
steady-state,  laminar  solution  for  a  given  flow  geometry  and 
a   small   perturbation   flow   which   is  added  to  the  laminar 

flow.   The  velocity  and  vorticity  become,  respectively,   V+v 

and  ft+aj   where  v  and  u;  are  the  perturbation  guantities.   It 

is  easy  to  show  that  the   perturbation   flow   must   satisfy 
continuity  independently. 


vv   =  0 


It  can  also  be  shewn  that 


(A-7) 


u   =    vxv  . 


(A-8) 


Substituting   V+v   and  ft+w  into   the  vorticity   transport 
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equation  and  subtracting  the  equation   for   the   unperturbed 
flow,  the  result  is 


1_V2w   ♦    vxtVXtu-^Xv  +  vxoj)       -    dui    =    0    .  (A-9) 

"Re  "31 

This   equation   is   linear   in   the   perturbation  quantities 

except  for  the  second  order  term  vxw.   If   the  perturbation 

is  small,  this  term  can  be  neglected,  linearizing  the 
equation.  The  result  is  the  linearized  perturbation 
vorticity  transport  equation. 

1  v2a>  +  VX(VXw-(ftXV)   -  3a;   =0  (A-10) 

"Re  d^ 

A  second  form  of  this  equation  is  derived  by  expanding 
the  curl  cf  the  cross  product  terms  and  noting  that  v»V  = 
v«v   =    v£l      =   voj   =    0.      The    result    is 

(1/Re)  v2<^    +     (cu«v)  V   -{V«v)oj    +     (v*v)ft 

-     {&•*)  v    -    duj/at    =    0    .  (A-11) 

A  third  form  of  eguation  (A-10)  is  obtained  using  the  vector 
identity 

v2w  =  v(v*w)  -  vx(vxcj)  (A-12) 

and   noting   that  w   is   non-divergent. 

-1    VXVXflJ   +    vXVXcj  -    VX&Xv       -    Soj    =    0  (A-13) 

He  3^ 
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APPENDIX  B 

DEPENDENCE  BETWEEN 
COMPONENTS  OF  THE 
PERTURBATION  VORTICITY  TRANSPORT  EQUATION 

The  vorticity  transport  equation 

-1_vxv*u>  +  vxVxw  -  vx&xv   -  3oj  =  0  (A-13) 

E"e  3t 

appears  superficially  to  represent  a  set  of  three,  coupled, 
simultaneous  equations.  This  equation  actually  has  only  two 
independent  conditions.  To  show  this,  equation  (A-13)  is 
first  simplified  by  the  introduction  of  the  velocity   vector 

potential  W  defined  by 

v  =  vxW  .  (2-7) 

In  Appendix  D  it  is  shown  that  this  results  in  a  vector 
equation  which  contains  only  the  three  components  of  W  as 
its  unknowns. 

-1_vxvxvxvxW  +  vxVxvxvxW  -  vxfixvxW 
Ee 

-  vxvxj>w  =  0  (D-1) 

Furthermore,  the  followinq  form  of  solution  is  assumed. 

X 
W(x,y,z,t)  =  £if(y)  +  jq(Y)  +  kh  (y)  ]e  (2-9a) 

where 

X  =  ax  +  /3z  +  yt      .  (2-10) 

Consider  equations  (D-1)  expressed  in  abbreviated  form  by 
the  statement 
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r  = 


*v 


V  *zJ 


=    0 


(B-1) 


which  may  be  expressed  as  three  separate  equations, 


r* 

0 

ry 

= 

0 

r 

z 

= 

0 

(B-2a) 


(B-2b) 


(B-2c) 


Equations   (B-2)  are  written  for  cartesian  coordinates.   The 

components  T  ,     T     and  T     happen  to  be  the   three  components 
x    y      z 

of   a   vector   obtained  by  applyinq  the  curl  operator  to  the 

momentum  transport  equation.   It  is  a  vector   identity   that 
the  curl  cf  any  vector  is  nondivergent  and  thus 


iL(T)   +^_(D  +  a_{D=o 
dx   x     3y   y    3z   z 


(B-3) 


With  the  assumed  form  of  W  given  by  equations  (2-9a)  and 
(2-10)  and  the  fact  that  all  these  equations  are  linear  in 
H#  (B-3)  becomes 

ar    +  i-(T)  +  &r  =0  (b-4) 

x   3y   y      z 

Equation  (B-4)  has  been  verified  for  the  resulting  form  of 
the  eguation  developed  in  Appendix  D  (equations  (D-46) 
through  (D-55) )  by  differentiating  the  second  eguation, 
multiplying  the  first  by  a  and  the  third  by  /3,  and  adding 
the  three  equations  thus  obtained.  All  terras  are  found  to 
cancel  out.   From  equation  (B-*l)  it  may  be  inferred  that   it 
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is  sufficient  tc  satisfy  the  second  of  equations  (B-2)  and 
any  linear  combination  of  the  first  and  third.  If  the 
second  is  satisfied,  then 


r    =  0 
y 

and  it  must  be  true  that 
|_(T)  =  0  . 

Eguaticn  (B-4)  reduces  to 

ar    +  pt   =  0  . 

x     z 


(B-5) 


(B-6) 


(B-7) 


As   discussed   in   section   2,  the  proper  choice  of  a  second 
condition  is 


■pr    +  aT  =  0 
x     z 


(2-11) 


(B-8) 


Equations  (B-7)  and  (2-11)  may  be  expressed  in   matrix   form 
as 

"  a  (3 
■j3     a 

If  the  determinant  is  non-vanishing,  that  isr  if 

a2  +  /32  *    0  , 
then  the  only  solution  of  equation  (B-8)  is 


(B-9) 


(B-10) 


Thus,  equations  (B-5)  and  (B-10)  verify  that  equation  (B-1) 
will  be  satisfied  by  the  two  equations  (B-5)  and  (2-11). 
That  is,  these  two  conditions  are  sufficient. 

For  pipe  flow  in  cylindrical  coordinates,  equation  (D-1) 
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may  be  expressed  similarly  to  (B-1)  as 


r  = 


>  =  0 


e ) 


which  are  three  separate  equations. 


r  =  o 

X 


(B-11) 


(B-12a) 


r   =  o 

r 


(B-12b) 


r   =  o 

e 


2-  ( r )  +  ji>_  (r  r )  +  n^iD  =  o 

ax      x    rar   r    ra5  6 


(B-12c) 


As   in    the   cartesian   case,   equation   (B-11)   must   be 
non divergent. 


(B-13) 


With  the  assumed  form  of  W   given   in   equations   (E-1)   and 
(E-2)  ,  eguaticn  (B-13)  becomes 

ar  +  jd (r  r  )  +  gr  =  o  (b-iu) 

x   r    r    re 

Eguation    (E-14)    has   been   verified   by   performing   the 

indicated  operations  on  T  ,     T  ,    and  T   .    All   terms  cancel 

x   r      e 

out. 


From  equation  (E-14)  it   may   be   inferred   that  it   is 

sufficient   to   satisfy   equation   (B-12b)   and   any  linear 

combination  of  T     and  T  .   If  the  second  is  satisfied  then 

X        G 


r   =  o 
y 


(B-15) 
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Therefore 


1D(r  D  =  0 
r    r 


(B-16) 


for  any  value  of  r  except  at  r=0,  where  it  is   true   in   the 
limit  as  r  approaches  zero.   Equation  (B-14)  then  reduces  to 


ar  ♦  £r  =  o 

x    re 


(B-17) 


As   discussed   in   section  2,    the  proper  choice  of  a  second 
condition  is 

-§T    +  aT  =0  (2-14) 

r  x     9 

These  two  eguations  may  be  expressed  in  matrix  form  as 


a   £ 

r 

~i      a   Jr  \         I  0\  (B-18) 

r 


and,   as  before,  equation  (B-11)  will  be  satisfied  under  the 
condition  that 


°z  +fl 


*    0 


(B-19) 
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APPENDIX  C 

REDUNDANCY  OF  THE  THREE  COMPONENTS 

OF  THE  VELOCITY  VECTOR  POTENTIAL 

AND  THEIE  RELATIONSHIP  TO  THE  STREAMFUNCTION 

1.   Cartesian  Coordinates 


To   find   the  relationship  between  the  components  of 
the  velocity   vector   potential   in   cartesian   coordinates, 


consider  the  definition  of  W. 


<^_    9_         d_ 
5x   ay   ^z 

f(y)  9(y)  h(y) 
which  has  the  components 
u(y)  =  Dh  -  /3g 


i  j  k 
a  D  yS 
f  g   h 


(c-1) 


(C-2a) 


v(y)  =  /9f  -  ah 


(C-2b) 


w  (y)  =  ag  -  Df  . 


(C-2c) 


The  components  of  the  velocity,  u,  v,  and  w,  are  not 
independent  hut  are  related  by  the  continuity  equation,  that 
is,  by  the  condition  of  non-divergence. 


u  =  au  +  Dv  +  /3w  =  0 


(C-3) 


Let  f (y)  be  identically  zero.   Then 
u(y)  =  Dh  -  £g 

v(y)  =  -ah 


(C-4a) 
(C-4b) 
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«(y)  =  <*g  -  (c-4c) 

Thus  the  components  of  the  velocity  may  be  represented  by 
the  two  components,  h  and  g,  of  the  velocity  vector 
potential  without  any  loss  of  generality  other  than  that 
imposed  by  the  condition  of  non-divergence.  If  g (y)  is  set 
to  zero,  the  following  set  of  eguations  result. 

u(y)  =  Dh  (C-5a) 

v(y)  =  £f  -  ah  (C-5b) 

w(y)  =  -Df  (C-5c) 

And  if  h  (y)  is  identically  zero 

u(y)  =  -/3g  (C-6a) 

v(y)  =  £f  (C-6b) 

w(y)  =  ag  -  Df  .  (C-6c) 

In  the  classical  case  of  plane  perturbations  the 
parameter  yS  is  zero.  Then  the  3-dimensional  eguations, 
(C-2)  ,  reduce  to 

u(y)  =  Dh  (C-7a) 

v(y)  =  -ah  (C-7b) 

w  (y)  =  ag  -  Df  .  (C-7c) 

From  eguations  (C-7)  it  can  be  seen  that,  for  the  case  of 
£=0,  letting  h(y)=0  cannot  be  allowed  because  both  u(y)  and 
v (y)  are  forced  to  be  zero.  Since  it  is  desirable  to  have 
the  three-dimensional  representation  of  the  problem 
reducible  to  the  classical  two-dimensional  case,  it  is 
stipulated  that  the  function  to  be  eliminated  not  be  h(y). 
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Note  that  the  classical  case  of  two-dimensional 
plane  perturbations  corresponds,  in  the  present  notation,  to 
setting  f (y)  and  g  (y)  egual  to  zero  and  retaining  h  (y) 
alone.  This  component  of  the  velocity  vector  potential 
corresponds  to  the  classical  definition  of  the  stream 
function  for  two-dimensional  flow.   From  equations  (C-7) 


u(y)  =  Dh 


(C-8a) 


v  (y)  =  -ah 


(C-8b) 


w(y)  =0  (C-8c) 

for  the  classical  two-dimensional  flow. 
2.   Cylindrical  Coordinates 

The  relationship  between  velocity  and  the  vector 
potential  in  cylindrical  coordinates  is  similar  to  that 
relationship  just  developed  for  cartesion  coordinates.    The 

definition  of  W  in  cylindrical  coordinates  is 


1e 
r   x 

1e 

r  r 

e 
e 

=    VXW    = 

a 

D 

£ 

f 

g 

rh 

which  has  the  components 

u(r)  =  JD(rh)  -  gg 
r        r 


(C-9) 


(C-10a) 


v(r)  =  it   -  ah 

r 


(C-10b) 


w(r)  =  ag  -  Df  (C-10c) 

As   before,   the   components  of  velocity  are  related  by  the 

continuity  equation 
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v«u  =  au  +  ID  (rv)  +  £v  =  0  .  (C-11) 

r        r 

It  is  easily  shown  that  the  velocity  components  in 
cylindrical  coordinates  also  can  be  represented  by  two 
components  of  the  velocity  vector  potential,  without  loss  of 
generality  other  than  that  imposed  by  the  condition  of 
nondivergence.  In  the  case  of  /3=0  and  h  (r)  =  0,  the 
eguations  become 

u(r)  =0  (C-12a) 

v(r)  =0  (C-12b) 

w(r)  =  ag  -Df  (C-12c) 

which  cannot  be  allowed.  Boundary  conditions  in  cylindrical 
coordinates  impose  the  following  restriction  on  P. 

/3  =  in,      n  =  0,1,2,...  (C-13) 

The  case  of  n=0  is  considered  separately  from  n*0.  For  n=0 
it  is  necessary  to  retain  h  (r)  and  set  either  f (r)  or  g (r) 
egual  to  zero.  For  the  other  cases,  n*0,  h  (r)  may  be  the 
component  of  the  vector  potential  which  is  set  to  zero. 
Note  that  n=0  corresponds  to  axially  symmetric 
perturbations.   For  f  (r)  =0 

u(r)  =  ID(rh)  -  @>g  (C-14a) 

r        r 

v(r)  =  -ah  (C-14b) 

w(r)  =  ag  .  (C-14c) 

For  g (r)  =0 

u(rj  =  JD(rh)  (C-15a) 

r 
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v(r)    =  if   -ah 

r 


(C-15b) 


v (r)    =   -Df 
And   for    h  (r)    =0 

u(r)    =   -£g 

r 


(C-15c) 


(C-16a) 


v(r)    =  £f 

r 


(C-16b) 


w  (r)    =   ag    -Df 


(C-16c) 
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APPENDIX    D 

DEEIVATION    CF    COEFFICIENTS    FOR 
PLANE    POISEUILLE    FLOW 


The   basic   eguations   are 


-1_vxvx<*>  +  vXVxw-vxftxv    -   dui/dt   =    0  (A-13) 

I?e 


w   =    vxv  (A-8) 


V    =    i3(1-y2)    =    iu(y)  (2-3) 


&  =    k3y  (2-4) 

v   =    vxW  '  (2-7) 

The   vector   functions   W,   v   and  a)   are  assumed   to  have 

solutions  of  the  form  of  eguations  (2-9) . 

-      -      -      X 
W(x,y,z,t)=  [if  (y)+jg(y)+*h(y)  je  (2-9a) 


X 
v(x,y,z,t)  =  [iu(y)+jv(y)+kw(y)  ]e  (2-9b) 


-     -      -      X 
«(x,y,z,t)  =  [i£(y)+j^(y)+k£(y)  ]e  (2-9c) 

where 

X  ~  ax    ♦  /8z  +  yt      .  (2-10) 

Equation  (A-13)  can   be   written   in   terms   of   only   W   by 
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substituting  for  o>  using  equation  (A-8)  and  then 
substituting  for  v  using  equation  (2-7) .   The  result  is 

-1_VXVXVXVXW     +     VXVXVXVXW    -     tfX^XVXW 

"Re 

-  VXVX^W  =  0  (D-1) 

From   the   assumed   form   of   the   solution  for  the  velocity 

vector  potential  W,  given  in  equations  (2-9a)  and  (2-10),  it 
can  be  seen  that  the  partial  derivatives  d/dx,  a/3z,  and 
d/3t  become  multiplication  by  a ,  /3  and  y,  respectively.  The 
partial  derivatives  d/dy,  92/t=)y/  etc.  will  be  represented  by 
the  operator  symbols  D,  D2,  etc. 

The   easiest   method   for   determining  the   form  of  the 
coefficients  of  eguation  (D-1)  is  to  represent  the  operators 

as   matrices   multiplying   the  vector  W  and  its  derivatives. 

Brackets  are  used  to  indicate  a  3  by  3  matrix.  Using  braces 
to  enclose  the  components  of  a  vector,  the  following 
definitions  are  assumed; 


(D-2) 


(D-3) 


Higher  order  derivatives  are  defined  similarly. 

It  is  necessary  to  determine  the  matrix  operator  [ CP  ] 
which  can  be  used  to  replace  the  cross  product  operator, 
that  is,  a  matrix  [CP]  which  satisfies  the  relation 

[CP]K  =  tx[s]W  (D-4) 
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where   t    is   an   arbitrary      vector      and      [s]      is      an      arbitrary 
matrix.         If      the      components   of    t   are   t    .    t      and   t      and  the 

12         3 


components  cf  [s]  are  s  ,    then 

ij 


'** 


[CP]W  =  U\x 


t 

\     3/ 


s   s   s 

11    12    13 


s   s   s 

21    22    23 


s   s   s 

3  1    32    3  3 


(D-5) 


The  matrix  [CP]  may  be  obtained  by  multiplying  vector  W  by 
matrix  [s],  then  finding  the  cross  product  of  the  resulting 
vector  with  vector  t  and,  finally,  collecting  terms  in  f,  g, 
and  h  to  separate  out  the  vector  W.   The  result  is 


[CP]  = 


(t  s   -t  s   )  (t  s   -t  s   )  (t  s   -t  s   ) 

2  31    3  21  2   32    3  22  2  33    3   23 

(t  S    ~t  S    )  (t  S    ~t  S    )  (t  S    -t  S    ) 

3  11    1   31  3   12    1   32  3  13    1   33 

(t  S    ~t  S    )  (t  S    -t  S    )  (t  S    -t  S    ) 

1   21    2  11  1   22    2   12  1   23    2   13 


(D-6) 


A  matrix  operator  is  also  needed   to  replace   the   curl 
operator,  vx. 


f  (y)  X 

vxH    =    vx}<g  jyji  e 


i        J        * 

«L     «L     ^_ 

5x     3y     "3z 

XXX 
fe      ge      he 


—  X        -  X        -  X 

=   i(Eh-/%)e      +    j(/3f-cth)e      +    k(org-Df)e 


(D-7) 
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This  can  fce  written 


vxW    =    [A]W    +    [B]DH 


(D-8) 


where  [A]  and  [B]  can  be  determined  by  collecting  the  proper 
terms. 


[A]  = 


0 

-£ 

o" 

/3 

0 

-a 

0 

a 

0_ 

(D-9) 


[B]  = 


0 

0 

r 

0 

0 

0 

-1 

0 

0 

(D-10) 


This  operation  defined  by  equation  (D-8)  holds  for  the   curl 

/  \  x 

of  any  vector  with  the  form  (\i  (y)  ,v  (y)  ,w  (y) )  e  ,    that  is,  for 

any  curl  operation  performed  in  equation  (D-1). 

To   evaluate   the   third  terra  in  equation  (D-1)  take  the 

cross  product  of  -Q,  and  vxW,  that  is, 

JO  >X(£AJW  *  [BjDW)  . 

This  can  te  evaluated  using  eguations  (D-6)  and  (D-5) .    The 
result  is 


ftxvxW  =  [C  ]H  +  £C  ]DW 

1  2 


(D-11) 


where 


[C,] 


-3/3y    0   3ay 
0   -3)ey   0 

0     0    0 


(D-12) 
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I.C,]- 


0 

0 

0  " 

0 

0 

3y 

0 

0 

0 

(D-13) 


Next   operate    with   -vx 


•VX([C    ]W  +  £C    ])DW 


~    CA](£Ci]W+[C2]DW)     -    [B]D([Ci]W+i;C2]DW) 


=   "   [A]£C     ]W   -    £AJ£C    JDW    -    £B]£C     ]DW 


-  2-  _ 

-    [B]E£C     ]W    -    [B][C     ]D    W    -    £B]B[C     ]DW 


(    "    [AjECJ    «    [B]D£Ci])W 


+     ("£A][C2]    -    [BJtCj    -    £B]D£C2]DW 


2- 

"    £E]£C2]D    H 


(D-14) 


Note  that  B£ C  ]  means  9   £C  ].   The   coefficient   of   D2W   is 
1        3y   l 

egual  to  the  zero  matrix,  so  (D-14)  may  be  written 


-VX&xvW  =  [E  JW  +  £E  ]DW 

1  2 


(D-15) 


The   matrices   [E  ]   and  £E  ]  may  be  found  by  simple  matrix 
multiplication  and  addition. 


£B    ]   = 


0      -3&*y         0 
3/32  y        0        -3a/3y 

-3/8     3a/3y      3a 


(D-16) 


9a 


•£»,]- 


0      3£y 
0        0 


-3/3y      0        0 


(D-17) 


To      calculate      the   other   terms,    it   is   necessary    to    first 
find   vxvxW 


vx  (vxW)    =    vx(£A]W   +   [B]DW) 

=[A][A]W    +    [A][B]DW    +    [B][A]DW 

+    [B]C£A]W    +    [B][B]D2W    +    [  B  ]D[  B  ]DW 
This    may    he    written 


(D-18) 


VXVXH    =    [E     ]W    +    [E     ]DW    +    [E     ]D2W 

3  4  5 


(D-19) 


where  the  matrices  [E  ],  [E  ]  and  [E  ]  are   determined   from 

3*  5 

the  coefficient  matrices  of  equation  (D-18) . 


[E3]  = 


-/32     0     a/3 
0   -(a2+£2)   0 


a/3 


0 


-a2 


(D-20) 


[«J 


0 

a      0" 

a 

0     /3 

0 

$    o_ 

(D-21) 


[E  ]  = 

5 


-1 

0 

o" 

0 

0 

0 

m  0 

0 

-1 

(D-22) 


The  term  -vxvx^M  can  now  be  evaluted. 
31 
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■vxvx<5w   = 

31 


-yVXVXU 


=    -y([E3]W    *    [EJDW    ♦    [Es]D2W) 


(D-23) 


To  evaluate  the  term  vxVxvxvxW,  first  find  the  matrix 
operator  that  is  eguivalent  to  the  cross  product  of  V  and 
vxvXW.   Using  eguations  (D-19)  and  (D-6) 


Vx (vxvxW) 


-\rh 


E     ]W+[E     ]DW+[E    ]D2W) 

3  c      4  5 


=  rc   ]w  +  re   ]dw  +  re   id2w 

3  ♦  5 


(D-24) 


where 


<V 


0 
-a/3U 


0 

a2u 


0    -  (a2+/32)  u    0  _ 


(D-25) 


ICJ- 


0       0         0 

o  -/3u     0 

_aU       0      £U_ 


(D-26) 


£Cs]  = 


Now  operate  with  vx. 


0 

0 

o" 

0 

0 

u 

p 

0 

0. 

(D-27) 
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VXVXVXVXW    =    VX(£C     lW-t-fC     ]DW4[C     ]D2W) 

3  ♦  S 


(£A]+£B]D)  ([C     ]W+[C     JDW+£C     ]D2W) 

3  4  5 


=    £E     ]W    ♦    £E     ]DH    +    [E     JD2W    +    [E     ]D3W  (D-28) 

6  7  8  9 


where 


o 


a/32U    3  (a2+^2)  y   -a2/3u 
0         a(a2+/92)  u        0 


-a2#U 


a3U 


(D-29) 


[*,] 


"-3  ay    -a2u   -3/3y 
-a2U        0      -a/3U 

.  0      -a£u         0 


(D-30) 


tE8]  = 


"aU 

0 

o" 

0 

0 

0 

.  0 

0 

aU_ 

(D-31) 


[E    ]  =   £0],    the   zero    matrix. 

9 


(D-32) 


To  evaluate  the  term  -1_vxvxvxvxW,  first  apply  the   curl 

Ee 

operator  to  vxvxw  which  is  given  by  equations  (D-19)  through 

(D-22)  . 


VX(VXVXW)     =     (£  A]+[B]D)    ([E     ]W+£E     ]DW+£  E    ]D2W) 

3  4  5 


=    £C     ]W    +    £C     ]DW    +    £C     ]D2W    +    [C     ]D3H  (D-33) 

•■6  7  8  9 
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where 


£V  = 


/3{Cf2+/32) 

0 


0 
-  a(a2+/?2; 


0 
a(a2+/32) 

0 


(D-34) 


[C  ]  = 

7 


0      0   -(a2+yS2) 
0      0       0 


(a2+/32)   0 


0 


(D-35) 


tCe^ 


'  0 

£ 

o" 

-£ 

0 

a 

0 

-a 

<L 

(D-36) 


0 

0 

-1" 

0 

0 

0 

1 

0 

0 

•4 

ic,]  = 


Finally,  operate  once  more  with  the  curl. 


(D-37) 


VX(VXVXVXW)     =     ([A]+[  B]D)    ([C     ]W+[C     ]DW+[C     ]D2W+[C     jD^W) 

6  7  8  9 


=    [E        ]W    +    [E        ]DW    +    [E        1D2W    +    [E       ]D3W 

1      10  11  12  13 


♦    [E       ]D*W 

14 


(D-38) 


where 


[E       ]    = 

10 


/32(a2+/32j  o  -a/3(a2+/32)' 

0  (a2+y92)2  o 


-a£(a2+£2) 


a2  (a2+yS2) 


(D-39) 
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£E       ]    = 
11 


[E       ]    = 

12  J 


0  -a(a2+yS2)  0 

-a(a2+£2)  0  -/3(a2+/32) 

0  -/3(a2+/32)  o 


a2  +  2yQ2  0  -a/3 

0         (a2+£2)  0 


-a/3  0  2a2+£2 


[E       j    = 

13 


0 

-a 

o" 

-a 

0 

-£ 

0 

-£ 

0 

[E       ]    = 

14 


The  final  eguation  becomes 


"1 

0 

o" 

0 

0 

0 

_o 

0 

1_ 

-1    [E       ]D*W    -    1    [E       1D3W    +     £-1    [E       ]+[E    ])D2W 

+     (-1    [E       ]+[E     ]+[E     Y)  DW 
V  ^gL     iiJL     7JL     2J/ 

+     (-1    £E       ]+[E    ]+[E    ])W 
N  Ee      10  6  l  ' 

~    /([E     ]D2W+[E     ]DW+[E     ]W)     =    0 


which  can  be  written 


[M     ]D*H    +    £M     ]D3W    +    [M     ]D2W    +    [M     ]DW    +    [M     ]W 
+    yft  N     ]D2W    +    [N     ]DW    +    [N     ]w)     =0 


(D-40) 


(D-41) 


(D-42) 


(D-43) 


(D-44) 


(D-45) 


where 
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[MJ 


-10       0 

He 

0       0      0 


0   -1 

He 


(D-46) 


[M3] 


0       a    0 


<L    0       8 

Ee  Re 


(D-U7) 


[M2J 


-/32  +  T 
He 

0 


He 


-1       (0f2+£2) 

He 


He 
0 


•a2+T 
He 


(D-48) 


£Mi3 


-3ay    -aT      0 
-aT         0      -/ST 

-3/?y    -£T      0 


(D-U9) 


[H,] 


£2T  3a2y  -a/5T 

3/32y  (a2+#2)T      -3a£y 

-a/Sl-3/3  3a/3y  3a+a2T 


(D-50) 


[NJ 


1 

0 

o" 

0 

0 

0 

0 

0 

1 

(D-51) 


[NJ 


0 

-a 

o" 

a 

0 

-P 

0 

-/3 

0_ 

(D-52) 
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Ci03- 


$z  0    -a£ 

0   (a2+/32)   0 


-ap         0 
and  T  is  defined  by 


T  =  aU  -  1_  (a2+£2)  = 
le 


3a(1-y2) 


1   (<j2+£2) 

Eg 


(D-53) 


(D-54) 


Finally,  equation  (D-45)  can  be  divided  through  by  e  to 
yield  the  results  below. 
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APPENDIX    E 

DERIVATION    CF    COEFFICIENTS    FOR 
PIPE    FLOW 


The   basic  eguations   are 


-1    vxvxvxvxw    +    vxVxvxvxW    -   VX^XVXW-  vxvx  3W   =0  (D-l) 

Ee  9t 

V   =    e   2(1-r2)    =   e   0  (r)  (2-5) 

xx 


&  =    e   4r  (2-6) 

e 


and   the   velocity   vector  potential  W  has  a  form  similar  to 

that  for  the  plane  flow  case. 

X 
W(x,r,e,t)  =  [e  f(r)+e  g(r)+e  h{r)]e  (E-1) 

x      r      9 

where 

X  =  ax  +  fie    +  yt    .  (E-2) 

As  in  the  plane  flow  case,  the  partial  derivatives  3/<3x, 
2/de  and  3/at  become  multiplication  by  a,  ft  and  /, 
respectively.  The  partial  derivatives  ^/^r,  dz/drz ,  etc. 
will  be  represented  by  the  operator  symbols  D,  D2,  etc. 

The  nethod  cf  solution  for  cylindrical  coordinates 
proceeds  exactly  as  that  for  the  cartesian  coordinate  case 
presented  in  Appendix  D.  Definitions  similar  to  (D-2)  and 
(D-3)  are  assumed. 

^(girlie*  (E-3, 
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DW  ,<|| 


(E-4) 


Higher  crder  derivatives  are  similarly  defined.  The  matrix 
[CP]  defined  by  eguation  (D-6)  is  valid  for  cylindrical 
coordinates  and  is  used  to  evaluate  cross  products.  The 
curl  operation  in  cylindrical  coordinates  is  found  below. 


«w  =  „  jjj  [| 


e  = 


1e 
r   x 

1e 

r   r 

e 
e 

d 

3x 

9 
Sr 

33 

X 
fe 

X 

ge 

X 
rhe 

—  Y      —  Y     —  Y 

=  Je  (h+rDh-/3g)e   +  Je  fef-rah)e   +  e  (ag-Df)e    (E-5) 


r  x 


This  can    be    written 


r   r 


VXW    =    £A]W    +    [B]DW 


where 


£A]   = 


r 


& 

r 


0      -a 


(E-6) 


(E-7) 


[B]    = 


0 

0 

1~ 

0 

0 

0 

1 

0 

0_ 

(E-8) 


Following    the    method   of    Appendix    D,    the   first   step    is  to 
use    eguaticn     (D-6)     to   evaluate  ftxvxW. 
'0  ' 


io  Vx([A]+[b]dk)      =  [cjw  +  [cz]d;; 


(E-9) 
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where 


[C,]   = 


-4/3      0      4ar 
0      -4/3     4 

0         0         0 


(E-10) 


IC,]- 


0 

0 

0 

0 

0 

4r 

0 

0 

0_ 

(E-11) 


Next,    operate    with   -vx 


-vx&xvxW    =    -(£A]+£B]D)  ([C    ]W+[C    ]DW) 


=  £E    ]W    +    £E    ]DW    +    £  0]D2W 


(E-12) 


where 


0       -4/32 

r 

Hi 

r 

IEJ    = 

4/32        0      -4a/3 

r 

0        4a/3        0_ 

0      0      4/3 

£«,]- 

0      0      0 
-4£   0      0 

(E-13) 


(E-14) 


Next,    find    vxvxW 


VX(VXH)     =     (£A]+[B]D)   (£A]W+£B]DW) 


=    £E     ]W    +    £E    ]DW    +    £E     ]D2W 

L      3  ♦  5 


(E-15) 


where 
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IB3]  = 


0   - 


r 


a 

r 

<a2+g£) 
r* 


■A 


[BJ- 


a 
0 

r 


r 
r 


-1 

0 

0 

0 

0 

0 

0 

0 

-1 

— 

r 

-<X2  +  1 


I>.]- 


The  term  -vxvx<9w  can  now  be  evaluated. 
3^ 


(E-16) 


(E-17) 


(E-18) 


■vxvxSw    =    -VVXVKW 
9^ 


=    -r([E3]H    +    [EJDW    +    [E5]D2WJ 


(E-19) 


Next  find  VxvxvxW  using  equations  (E-15)  and  (D-6) 
Vx(vxvxW)  = 


j  U(jr>(  X  ([  E3  ]i*[  E^  ]DW+[  Es  ]D*W) 


=  re   ]w  +  re   ]dw  +  re   ]D2w 

3  4  5 


where 


- 
0 

0 

- 
0 

^^' 

-agU 
r 

r*~ 

(a2-1    )U 
r* 

0 

"(a2+£2)U 
r? 

4U 

r* 

(E-20) 


(E-21) 
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[c4] 


0       0         0 

0    -£U    JU 
r      r 


aU      0 


r 


(E-22) 


[C.J- 


0 

0 

— 

0 

0 

0 

0 

0 

0 

0 

(E-23) 


Now,    operating    with    vx,    the   result   is 


VXVXVXVXW    =    OX([C     ]W+[C     ]DW  +  [C     ]D2W) 


=   [E     ]W    +    £E    ]DW    +    [E     ]D2W    +    [E     lD^W 

6  7  8  9 


(E-24) 


where 


a£2u      4r  (a2+^2)  -a£D   -4g2-ga2U 
r?  r7     r  r        r"~ 


0 

-a£^D 

r 


a(az*g2)0 
r? 

ct/?U 
r* 


-aBU 

(a3-a)U 
r 


(E-25) 


CV   = 


aU-Ura  -a2u   -4£ 

r 


-a2D 
0 


0      -a£U 

r 


r 


aU 

r 


(E-26) 


aU 

0 

0 

0 

0 

0 

0 

0 

aU 

£E    ]« 


and      the      matrix     £E    ]      is      identically   zero. 

9 


(E-27) 

For  the  final 


terra,  first  calculate  vxvxvxW. 
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VX(VXVXW)     =     (£A]+£B]D)    ([  E     ]W+[E     ]DW+£E     ]D2W) 

<S  ^  9 


=    £C     ]W    +    £C     ]DW    +    [C     ]D2W    +    [C     ]D3W  (E-28) 


where 


£V 


r  r? 

-?£2 

r^ 


4+£(a2+££) 


2a 


a-a(a2+^2) 


"J   ~JL(a2+^2) 
r"3    r  r? 

-  a  +  a(a2+fi2) 
r*  r7 

r^ 


(E-29) 


0 

"4 

r^ 

[c7]  = 

"ft 

0 

(a* 

-1 
r* 

-a 

r 

0 

r 

-2 

r 

C<V   = 

r 

0 

a 

1 

-a 

0 

r 

-  I  nZ 


(a 


r7     r< 


a 

r 


(E-30) 


(E-31) 


*V    = 


And/  finally/ 


0 

0 

-1 

0 

0 

0 

1 

0 

0 

(E-32) 


VX(VXVXTXri)     =     (£A]+£B]D)    ([C6]W+[C7-DW  +  [C8]D2W  +  [C9]D3W) 

=    £E        ]H    +    £E        ]DW    ♦    £E       ]D2W    +    [E13]D3W 


+    £E       ]D*W 


(E-33) 
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Letting 
t  = 


(c*2+£|) 


(E-34) 


the  coefficients  are 

f 5"  r*  f* 

-a£t  ig+3fit 

r  r*    r^ 


[E       ]   = 
10 


-a£-a£t 
r^    r 

-A-2a*g-J|t 
-3   -2a2-3^2+a2t 


(E-35) 


[E        ]    = 
11 


""-3^2+1 
r"3    r^ 

ft 

r 

a  -at 

f"2 

-at  +  o 
f"2 

-£2+a2 

r7   r 

-a£ 
r* 

-3§-fit 

P   r 

r- 


r     r7 

3    +a£-2^2+t 
r"7  r        r^   r 


(E-36) 


[E       ] 

12 


*fi£-1    +t 

r"z    r"z 

-a 

r 

-a$ 

r 


-2a 

r 


2/3 


r 
-2/3 

d2-3    +t 
r* 


(E-37) 


EE       ] 

13 


2 

r 

-a 


-a 


0      - 


-i 

r 

2 

r 


(E-38) 


[E14]   = 


1 

0 

0 

0 

0 

0 

0 

0 

1 

(E-39) 


The    final    equation,    with    all    the    terms    grouped      together 

X 

has   the   same   form  as  equation  (D-43)  and,  dividing  by  e  , 
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can  be  written  in  the  form  of  equation  (D-55) . 

iD3f~)  (D2g) 


♦   /d«: 


[■.3ft?st  *  CM3^glg\  +  <£H2>*N2»/g24' 


-    C£Ha3*xC«43lJg|}  -    O0]*rEN0])[g 


=o 


(D-55) 


where 


[MJ   = 


-10      0 
Re 

0       0      0 


0      0-1 
Re 


(E-40) 


[Ml   = 


~_  2_        q_        0 
rRe      Re 


a 

Re 


0  £ 

rRe 

fie      rRe 


(E-41) 


C«23 


T+        1    - 
r^Re 

_a_ 

rRe 

a/3 

r~Re 


2a_ 

rRe 

-t 
Re 

•_2/3 
r^fie 


rRe 

i 

r7Ke 

T+       3    - 
r^Re 

■as 

Re 

(E-42) 


[MJ 


1T-4ra+3££_-    1 
r*  r^Re  r^Re 

-aT-     a 
r~zRe 

_ajS_-4y9 
r"2Re 


-aT-    a 
f^Re 


-  aB 

r^Re 


r^Re   rRe  r      r^Re 

-£T+    3/3        JT+2fi2_-_3      __a£ 
r      if^Re      r      r^Re   r7Re   rRe 

(E-43) 
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r*      T*"Re 

r      r^fie 

-g£T 

r 


-aT+4ra2+     a 
r  f^He 

tT+    az_-_P2_ 

"if^Ee   r*le 

^  T  +  4a£-_3/3-_2gt 
r"2  rTEe   r?Re 


-g£T+_qfi 

r        r"*Ee 

-J?T-4a£+     £     +2az@ 

(a2-_l    )T+    3    _+_a2_+    2g£ 

r"2        iTHte   r^He  pRe 


(E-44) 


[>,]  - 


1 

0 

— 

0 

0 

0 

0 

0 

0 

1 

(E-45) 


[Hf]    = 


x 

r 
-a 


•a 

0 

r 


r 

r 


(E-46) 


[Nq]    = 


r"2  r  r 


-4 


0  t 

r        r^  r 


(E-47) 


and 


aO  -   t      =    2a{1-r2)    - 
Se 


i    <«2+£|> 

Ee  r? 


(E-48) 


and   t   is    defined    fcy    equation    (E-34) 
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APPENDIX  F 


BOUNDARY  CONDITIONS  AT  THE  WALL 


1.   Boundary  Conditions  for  Plane  Flow  at  y=±1 

The  boundary  condition  is  that  the  velocity  v  is 
zero  at  the  walls.  The  boundary  conditions  for  W  must  be 
determined.   Ey  the  definition  of  W, 


V  =  VXH  . 


(2-7) 


The  form  of  v  and  W  is  given  by  equations  (2-9a) ,  (2-9b)  and 
(2-10).   Equation  (2-7)  can  be  written 


u(y)  e 


v  (y)  e 


w  (y)e 


/     x\ 
f  (Y)  e 


vx(g(y)e  J> 
X 


h(y)e 


i  J  k 

9/3x         ^/<3y  2/dz 

XXX 
f(y)e     g(y)e     h(y)e 


i  3  * 

a  D  P 

f(y)   g(y)  My) 


-  X       -  x       -  x 

=  i(Dh-£g)e      +    j(/3f-ah)e      +    k(ag-Df)e 


(F-1) 


Equation    (F-1)    must    be   satisfied   by   each   component 
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separately.   Dividing  by  e  ,  the  results  are 
u(y)  =  Dh(y)  -  £g(y) 

v(y)  =  /3f  (y)  -  ah(y) 


(F-2a) 
(F-2b) 


w(y)  =  ag(y)  -  Df(y)  (F-2c) 

As  described  in  Appendix  C,  one  component  of   the  velocity 

vector  potential  W  may  be  arbitrarily  set  equal  to  zero, 

a.   Case  1 f  (y)  =  0 


At  y  =  ±1,  equations  (F-2)  become 
u(±1)  =  0  =  Dh  (±1)  -  £g(±1) 

v (±1)  =  0  =  -ah (±1) 

w(±1)  =  0  =  erg  (±1)  . 
These  equations  imply 
g(±1)  =  0 

h(±1)  =  0 

Dh(±1)  =  0  . 

b.   Case  2. ..   g  (y)  =  0 

At   y    =   ±1,    equations    (F-2)    become 
u(±1)     =    0    =    Dh(±1) 

v(±1)     =    0    =  /3f  (±1)     -  ah(±1) 

w(±1)     =    0    =    -Df  (±1)     . 


(F-3a) 
(F-3b) 
(F-3c) 

(F-4a) 
(F-Ub) 
(F-Uc) 


(F-5a) 
(F-5b) 
(F-5c) 
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These  equations  imply 

/3f(±1)  =  oh(±1)  (F-6a) 

Df(±1)  =  0  (F-6b) 

Dh{±1)  =  0  .  (F-6c) 

c.   Case  3-  . .   h  (y)  =0 

At  y  =  ±1,  equations  (F-2)  become 
u(±1)  =  -/3qi±-\)  (F-7a) 

v(±1)  =  pf  (±1)  (F-7b) 

w(±1)  =  org  (±1)  -Df(±1)  .  (F-7c) 

These  eguations  imply 

f(±1)  =0  (F-8a) 

Df(±1)  =  0  (F-8b) 

g(±1)  =  0  .  (F-8c) 

2.   Boundary  Conditions  for  Pipe  Flow  at  r  =  1 

The  boundary  condition  is  that  the  velocity  v  is 
zero  at  the  wall.  The  boundary  conditions  for  K  must  be 
determined.   Ey  the  definition  of  W, 

v  =  vxW  .  (2-7) 

The  form  of  v  and  W  is 
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v(x,r,e,t)     =    (e    u(r)+e   v  (r)  +e   w(r))e 

x  r  e 


(F-9) 


W(x,r,e,t)    =     (e   f(r)+e    g  (r)  +e   h  (r)  )  e 

x  r  9 


(F-10) 


where 


X   =    ax    +  $6    +  yt    .  (F-11) 

In   cylindrical    coordinates,    equation    (2-7)    can   be   written 


/  x 

u  (r)  e 


w  (r)  e 


/, 


XI  X 

v(r)e     >  =    vx/g(r)  e    \ 


f  (r)e 


h(r)e 


Je     Je        e 
r   x   r   r        9 


<9        9        d 
3x"      3f     50" 

XXX 
fe      ge      rhe 


Je     J[e      e 

r   x   r   r      9 


D      0 


f        g      rh 


=  Je    (p  (rh)-£g)e      +  Je    (fef-arh)  e        +   e    (ag-Df) 

(F-12) 

Equation       (1-12)       must      be      satisfied      by        each        component 

X 

separately.   Dividing  by  e  ,  the  results  are 


u(r)  =  JD(rh(r))  -  fig  (r) 
r  r 

=  Jh(r)  +  Dh(r)  -  £g(r) 
r  r 


(F-13a) 


v(r)  =  £f  (r)  -ah(r) 

r 


(F-13b) 
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w(r)  =  ag(r)  -Df  (r)  (F-13c) 

As  described  in  Appendix  C,      one   component  of   W   may  be 
arbitrarily  set  equal  to  zero. 

a.  Case  1...   f  (r)  =  0 

At  r=1,  eguations  (F-13)  become 
u(1)  =  0  =  h(1)  +  Dh{1)  -  >9g  (1)  (F-14a) 

v(1)  =  0  =  -ah  (1)  (F-14b) 

w(1)  =  0  =  ag  (1)  (F-14c) 

These  eguations  iaply 

g(1)  =  0  (F-15a) 

h(1)  =  0  (F-15b) 

Dh(1)  =  0  (F-15c) 

b.  Case  2.  . .   g (r)  =  0 

At  r= 1 ,    eguations  (F-13)  become 
u(1)  =  0  =  h(1)  +  Dh(1)  (F-16a) 

v(1)  =  0  =  0f  (1)  -ah(1)  (F-16b) 

w(1J  =  0  =  -Df  (1)  (F-16c) 

These  eguations  imply 

/5f  (1)  =  ah  (1)  (F-17a) 

Df(1)  =  0  (F-17b) 
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h(1)  =  -Dh(1) 

c.   Case  3...   h (r)  =  0 

At  r=1,  equations  (F-13)  become 
U(1)  =  0  =  -£g(1) 

v(1)  =  0  =  ySf  (1) 

w(1)  =  0  =  ag  (1)  -Df  (1) 
These  equations  imply 
g(1)  =  0 

f(D  =  0 
Df(1)  =  0 


(F-17c) 


(F-18a) 
(F-1 8b) 
(F-18c) 

(F-19a) 
(F-19b) 
(F-19c) 
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APPENDIX  G 

BOUNDARY  CONDITIONS 
FOR  PIPE  FLOW  AT  r  =  0. 
PERIODICITY   OF  6. 

1.   Boundary  condition  for  0 

The  condition  imposed  is  that   functions   be   single 
valued  at  8  =  0  and  Q  =  2ir.      The  assumed  form  of  the  solution  is 

X 
W(x,r,e,t)  =  [e  f(r)+e  g  (r)  +e  h(r)]e  (E-1) 

x      r      6 

where 

X  =  ax  +  £e  +  yt    .  (E-2) 

Thus  the  €  dependence  may  be  separated  as  the  factor 

pe 

e    =  cos  (-i/36)  +  isin(-i£e)  (G-1) 

This  terra  must  be  single  valued  at  6=0  and  Q=2ir.  Both  the 
real  and  imaginary  parts  must  match.  This  results  in  two 
eguations. 

cos  (-ip2ir)     =  1  (G-2a) 

sin  (-i£27r)  =  0  (G-2b) 

These  conditions  are  met  only  for  arguments  of  sine  and 
cosine  which  are  even  multiples  of  ir,    that  is 

(-i/?27r)  =  2n7r,    n=0/1,2/...  (G-3) 

or 

$   =  ni,    n=0,1,2,...  (G-4) 
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2.   Single-Valuedness  at  r=0 

It   is  reguired  that  a  vector  function  with  the  form 
of  eguaticn  (2-9)   be   single-valued   at   the   origin,   r=0. 

Consider  a  vector  V   which  has  the  form 

n 


V   =  I  e  u  (r)  +  e  v  (r)  +  e  w(r)  e 
x  n        r  n        Q  n 


(G-5) 


where 


X  =  ax  +  in8  +  yt 


(G-6) 


Figure  G-l     Vector  Function  Rotated  <f> 
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Consider  figure  G-1  where  point  P  is  rotated  an 
angle  4>  about  the  origin  to  become  point  P',  both  at  a 
distance  r  from  the  origin.   The  unit  coordinate  vectors  are 

e  ,    e  ,  and  e   at  P  and  e1 ,    e1,  and  e'  at   P*  .    Vectors   e 
x    r       9  x    r       6  x 

and   e1   extend   out  of  the  picture  normal  to  e   and  e  .   At 
x  re 

point  P' g    V   becomes 
n 

X« 

V1  =  le'u  (r)  +  e'v  (r)  +  e'w  (r)  ]  e  (G-7) 

n    I  x  n       r  n       9  n 

where 

X1  =  ox  +  in  (6+4>)  +  yt    .  (G-8) 

Then 

X»     X  ini 

e   =  e  e  (G-9) 

and 


e«  =  e  (G-10a) 

x     x 


e'    =   e   cosiji    +    e   sin<J>  (G-10b) 

r  r  G 


e1    =    e   cos<Ji    -    e   sin<|>  (G-10c) 

6  6  r 

Substituting    (G-9)    and    (G-10)       into       (G-7)       and      rearranging 
terms   results    in    eguation    (G— 11)- 

V    =      e   u    (r)    +    e    (v    (r)  cos<$>  -    v    (r)sin<()) 
n  xn  rvn  n  ' 

-     .  v       X    in<{)  (G-11) 

+    e     (v     (r)sin4)    +    w    (r)  cos4>)     e   e 
6  v  n  n  ' 
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It  is  required  that 

lira  V   =  lira  V  .  (G-12) 

r-*-0   n   r-»0   n 

If  equations  (G- 1 1)  and  (G-5)  are  set  equal  to  each  other  at 
r=0,  then  each  component  of  this  vector  equation  must 
separately  satisfy  the  equality.  The  result  is  the 
following  three  equations. 

in<J> 
u     (0)    =    u     (0)e  (G-13) 

n  n 


in<f> 
v(0)    =    [v    (O)cos(J)   -w    (0)sin<J)]e  (G-14) 

n  n 


incfc 
v    (0)    =    [v     (0)sin<J>    +    w    (0)cos<J)]e  (G-15) 

n  n  n 

Equation  (G-13)  can  be  written 

in(J) 
u  (0)[  1-e    ]  =  0  .  (G-16) 

n 

ini 
Por   n=0,    1-e      is    identically   zero,   so   u  (0)   is 

in<t> 
unrestricted.   For  n#0,  1-e     is,  in  general,  not  equal   to 


zero   so   u  (0)=0   for  n#0.   Equations  (G-14)  and  (G-15)  are 
n 

rewritten 


in<J>  in<t> 

v    (0)  (1-cos4>e        )    +    w     (0)  (sin<J>e        )    =    0  (G-17) 

n  n 


ind)  in<J> 

v     (0)  (-sin<J>e        )    +    w    (0)  (1-cos$e        )     =   0  (G-18) 

n  n 
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In  (G-17)  and  (G-18),  v  (0)  and  w  (0)  must  equal  zero  unless 

n         n 

the  determinant  of  the  coefficients  vanishes.    Setting   the 

determinant   of   the   coefficients   equal  to  zero  yields  the 

following  eguaticn. 

in<j>  in<J> 

(1-cos$e         )2    +     (sin$e        )2    =    0  (G-19) 

Multipling  the  terms  out,  using  Eulers  equation  and  noting 
that   sin2<J)  +  ccs2<j!=  1,    results   in 

1    -    2cos<J)  (co£n<J>+isinn<$)    +   cos2n<t>   +    isin2n<|>   =   0         (G-20) 

The  real  and  imaginary  parts  of  equation  (G-20)  must 
separately    he    satisfied   as   follows. 

1    -    2cos$cosn<J>    +   cos2n(J>    =    0  (G-21) 

-2cos<t>sinn$    +    sin2n<J>   =    0  (G-22) 

Using  well  known  identities  for  cos2n<J>  and  sin2n<|>  yields  the 
following    expressions. 

1    +   cosn4>  (co£n$-2cos<t>)    -sin2n<})   =    0  (G-23) 

sinn$  (cos<J)-cosn<j))    =   0  (G-24) 

For  n=0/  equation  (G-24)  is  satisfied  identically,  but 
equation  (G-23)  reduces  to 

2(1-cos$)  =  0  (G-25) 

which,  in  general,  is  not  satisfied  for  arbitrary  <t>.  For 
n=1,  equations  (G-23)  and  (G-24)  become 

1  -  (sin2(J>  +  ccs2<t>)  =  0  (G-26) 

sin<|!  (cos<J)-cos4))  =  0  (G-27) 

which  are  satisfied  identically  for  any  value  of  4>.  For 
n>1,  the  equations  will  not  be  satisfied.  In  particular, 
equation  (G-24)  requires  that  either  sinn<J)=0   or   cos<t>=cosn<J> 
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which  cannot  be  satisfied  for  arbitrary  4>  and  for  n>1. 


A  relationship  between  v  (0)  and  w  (0)  is  implied  by 

the   previous   result.    Using   equation   (G-17)  to  find  the 
ratio  of  v  (0)  and  w  (0)  results  in 

v    (0)  sm<J)e  sm<|> 


i<J)  -i<J> 

w    (0)  1-cos$e  e        -cos<J> 


-sin<J> 

=    l 


cos  (-$>)    +    isin  (-$)     -cos<J>  (G-28) 


or 


vt(0)    =   iw^O)  (G-29) 

Equation  (G-18)  yields  the  same  result.  A  summary  of  the 
results  of  applying  the  restriction  that  a  vector  function 
of  the  form  specified  by  equations  (G-5)  and  (G-6)  be  single 
valued  at  r=0  fcllcws. 

n=0   u  (0)  unrestricted  (G-30a) 


v  (0)  =  0  (G-30b) 

o 


w  (0)  =0  (G-30c) 


n=1   u  (0)  =  0  (G-31a) 


v  (0)  =  iw  (0)  (G-31b) 

l         i 
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n>1   u  (0)  =0  (G-32a) 

n 


v  (0)  =  0  (G-32b) 

n 


u  (0)  =  0  (G-32c) 

n 

3.   Symmetry  of  Solution 

Consider   figure   G-1   with  $>      =  tt.   From  equations 
(G-5)  and  (G-11) 

-  -  X 

V   =  [e  u  (r)  +  e  v  (r)  +  e  w  (r)  ]e  (G-33) 

n     x  n       r  n       9  n 


-     -         -         -        X   in7r 
V  =  [e  u  (r)  -  e  v  (r)  -  e  w  (r)  ]e   e  (G-34) 

n     x  n       r  n       0  n 

in-rr  n  - 

Note      that   e  =    (-1)    ,    and   let   the   components  of    V      and  V* 

n  n 

be  D  ,  V  ,  H   and  U',  V  and  W  •  .   From  eguations  (G-33)  and 
n    n    n      n    n      n 

(G-34)   the   components   of   V    and   V1   have  the  following 

n        n 

relations. 


U   =  U«  (-1)  (G-35) 

n    n 


n 
V   =  -V«  (-1)  (G-36) 

n     n 


W   =  -fi»  (-1)  (G-37) 

n     n 


Thus  the  components  of  a  vector  with  the  form  of  V   must  act 

n 
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either   as   an  even  or  odd  function  at  the  origin,  depending 

upon  n.   But   the   functional   dependence   of   U    on   r  is 

n 

contained   in   u  (r) .    Thus  u  (r)  must  look  like  an  even  or 

n  n 

odd  function.   Sinilarly,  v   and  w   must  behave  like  even  or 

n      n 

odd   functions.    In  summary,   eguations  (G-35) ,  (G-36)  and 

(G-37)  imply  the  following 


n  even     u  (r)    even  symmetry 
n 


(G-38a) 


'  (r)  \ 
n     >  o 


w  (r) 

n 


dd  symmetry 


(G-38b) 


(G-38c) 


n  odd     u  (r)    odd  symmetry 
n 


'  (r)  \ 
n     > 


w  (r) 
n 


even  symmetry 


(G-39a) 
(G-39b) 
(G-39c) 


4.   Summary 

If  the  curl  operation  is  performed  upon  a  vector 
with  the  conditions  of  symmetry  from  part  3,  the  resulting 
vector  obeys  the  same  rules  of  symmetry,  that  is,  the 
specification  that  even  or  odd  derivatives  of  the  components 
egual  zero  is  not  changed  by  applying  the  curl.  All  vector 
functions   are   assumed  to  be  of  the  form  of  equation  (G-7) . 

The  condition  of  single-valuedness  at  the  origin  for 
n=0  and  n=1  has  the  same  property.  Using  the  boundary 
conditions  at  the  origin  from  eguations  (G-30)  or  (G-31)  and 
applying   the   curl   operator   does   not  result  in  different 
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conditions  on  the  results.  This  can  be  verified  by 
performing  the  curl  operation  on  a  general  vector  of  the 
form  of  (G-7) ,  imposing  the  boundary  conditions  on  the 
original  vector's  components  and  observing  the  resulting 
restrictions  on  the  components  of  the  curl. 

For  n>1#  this  is  no  longer  true.  Only  the  cases  of 
n=0  and  n=1  are  considered  in  this  paper.  Thus,  for  n=0  and 
n=1,  the  vorticity,  velocity  and  velocity  vector  potential 
will  all  have  the  boundary  conditions  previously  discussed. 
The  boundary  conditions  on  the  velocity  vector  potential  at 
r=0  are  summarized  below. 

n=0   all  odd  derivatives  of  f  =  0 
all  even  derivatives  of  g  =  0 

(including  g  (0)  =  0) 
all  even  derivatives  of  h  =  0 

(including  h  (0)  =  0)  (G-40) 

n=1   all  even  derivatives  of  f  =  0 
(including  f  (0)  =  0) 
all  odd  derivatives  of  g  =  0 
all  odd  derivatives  of  h  =  0 
g  (0)  =  ih  (0)  (G-41) 
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c 

C  FFCGRAM  #1 

C 

C  PROGRAM  TO  PRINT  EIGENVALUES 

C  FOR  THE  3-D  POISEUILLE  FLOW  PROBLEM 

C 

c 

c 

c 

C      THIS  PROGRAM  SOLVES  THE  LINEARIZED  NAVIER  STOKES 

C      EQUATION  FOR  POISEUILLE  FLOW.   THE  EIGENVALUES 

C      RESULTING  FROM  THE  FINITE  DIFFERENCE  APPROXIMATION 

C      ARE  PRINTED. 

C 

C      INPUT 

C 

C  THE  FOLLOWING  MAY  BE  INPUT  TO  THE  PROGRAM  AS  DATA 

C  USING  NAMELIST,  'LIST'.   NOTE,  THE  DEFAULT  VALUES 

C  ARE  ONLY  SET  INITIALLY  AND  VALUES  SET  BY  THE 

C  USER  WILL  NOT  BE  CHANGED  BETWEEN  RUNS 

C 

C  N  -  HALF  OF  THE  NUMBER  OF  FINITE  DIFFERENCE  GRID 

C  PCINTS  ACROSS  THE  CHANNEL  NOT  INCLUDING  THE  END 

C  POINTS.   N  MUST  BE  . LE .  MDIM,  WHICH  IS  THE 

C  DIMENSION  OF  THE  MATRICES  IN  THIS  PROGRAM. 

C  DEFAULTED  TO  THE  VALUE  OF  NDIM,  THAT  IS  THE 

C  DIMENSION  OF  THE  MATRICES.   SEE  PROGRAM  BELOW  FOR 

C  THE  DEFAULT  VALUE. 

C 

C  REY    -    THE    REYNOLDS    NUMBER       (REAL*8)       DEFAULT 

C  VALUE    =    60C0.0 

C 

C  ARfAItBRfBI  -  THE  REAL  AND  IMAGINARY  PARTS  OF  THE 

C  PERTURBATION  WAVE  NUMBERS   <REAL*8) 

C  DEFAULTED  TO  O.O,  1.0,  0.0,  AND  0.0  RESPECTIVELY 

C 

C  VEL  -  THE  VELOCITY  OF  THE  MOVING  COORDINATE 

C  REFERENCE  SYSTEM  FOR  WHICH  THE  STABILITY  IS 

C  DETERMINED.   (REAL*3)   DEFAULTED  TO  0.0 

C 

C      CUTPUT 

C 

C  THE  CUTPUT  OF  THIS  PROGRAM  IS  A  TABULATION  OF  THE 

C  EIGENVALUES.   TWO  LISTS  ARE  PkINTED,  ONE  FOR  THE 

C  EIGENVALUES  CORRESPONDING  TO  EVEN  EI GENFUNCT I ONS 

C  ANC  CNE  FOR  THOSE  CORRESPONDING  TO  ODD  EIGEN- 

C  FUNCTIONS.   THE  STABILITY  OF  EACH  EIGENVALUE  IS 

C  PRINTED  AND  THE  LEAST  STABLE  EIGENVALUE  IS  MARKED 

C  WITH  ASTERISKS.   A  PLOT  OF  THE  EIGENVALUES  IS  ALSO 

C  PFINTED. 

C 

C      SUBROUTINES 

C 

C  THIS  PROGRAM  CALLS  THE  SUBROUTINE  'DEIGEO*  TO 

C  SOLVE  FOR  THE  EIGENVALUES.   SUBROUTINE  9PLOTP» 

C  IS  USED  TO  PLOT  THE  EIGENVALUES  ON  THE  PRINTER. 

C 

C 

C 

C 

IMPLICIT  REAL*8  (A-H,0-Z) 

CCMPLEX*16  A,B 

REAL*4  GR4(60) ,GI4(60) 

REAL*8  GRE(60),GIE(60) ,GR0(6  0) ,GI0(60) 

CGMPLEX*L6  XMAT(60,30,3) 

C0MPLEX*l6  YMAT(60,30) , WVECC60)  ,BMAT(5,60) 
EQUIVALENCES  MAT  (1,1)  ,XMAT(1,  1,3)  )  , 

*  (BMAT(  1,  1) ,XMAT(1 ,1,3)  ) , 

*  (WVEC(  1)  ,XMAT( 1,6,3)  ) 
NAMELIST  /  LIST  /  N  ,  RE  Y  ,  AR  ,  AI  ,  BR  ,  Li  I ,  VE  L 

C 

C  INITIALIZE  VARIABLES   (SET  DEFAULT  VALUES) 
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MDIM  =  60 

N  =  60 

REY  =  6CC0D0 

AR  =  0C0 

AI  =  ICC 

BR  =  OCO 

BI  =  OCO 
C 

C  READ  NAMELIST  AND  SET  ALPHA  AND  BETA 
C 

1  READ(5,LIST,END=1001 

A  =  DCMFLX(AR,AI) 

B  =  DCMPLX(BR,BI ) 
C 

C  PRINT  INPUT  VALUES  AS  PAGE  HEADING  FOR  EIGENVALUE  LIST 
C 

WRITE(6,9004) 

9004  FORMAT  ('  1«  ) 
WRITE(6,9005)  N, RE Y, A , B , VE L 

9005  FCRMATP  N  =',14,/,'  REY  =  ■  ,F10  .  2,  8X  ,  •  ALPH  A  =  ', 

*  2F12. 7f8Xi 'BETA  = •  ,  2F1 2. 7, / , •  VEL  =',F7.2) 
C 

C  CALL  SUBROUTINE  TO  SOLVE  FOR  EIGENVALUES. 
C 

CALL  DEIGEO(A,B,REY,N,MDIM,GRE,GIE,GRO,GIO,XMAT,YMAT, 

*  BMATtWVECj 
C 

C  DETERMINE  WHICH  EIGENVALUE  IS  THE  LEAST  STABLE. 
C 

TEMP  =  -1D10 

MARK  =  1 
C 

DO  20  1  =  1, N 

IF(GRC(I )+AR*VEL.LT.TEMP)  GO  TO  20 

TEMP  =  GRO(I )+AR*VEL 

ITEMF  =  I 
20  CONTINUE 
C 

DO  40  1  =  1, N 

IF(GRE( I )+AR*VEL.LT.TEMP)  GO  TO  40 

TEMP  =  GRE(I)+AR*VEL 

ITEMP  =  I 

MARK  =  2 
4C  CONTINUE 
C 

C  LIST  EIGENVALUES  FOR  ODD  E IGENFUNCTI ONS 
C 

WRITE(6,9003) 
9003  FORMAT*///, 6X, 'GAMMA  REAL1  , 5X, ■ GAMMA  I  MAG1  , 12X , • STA^' ) 

WRITE(6,9C06) 

9006  FORMAT! ^EIGENVALUES  FOR  ODD  EIGENVECTORS1,/) 
C 

CO  50  1=1, N 
TEMP  =  GPOU  )+AP*VEL 
WRITE (6, 9000)  GRC( I ) ,G 10 ( I ) , TEMP 
IF ( I .EG. ITEMP. AND. MARK. EQ.l)  WRITE (6,9001) 
5C  CONTINUE 
900  0  FORMAT! ■ 0*  , 1P2D1 5. 4, 1PD2  0. 4) 
9001  FCPMAT ( '  +  «  ,52X, •***• ) 
C 

C  LIST  EIGENVALUES  FOR  EVEN  EIGENVECTORS. 
C 

WRITE(6,9007) 
90C7  FCRMATCOEIGENVALUES  FOR  EVEN  EIGENVECTORS',/) 
DO  55  1=1, N 
TEMP  =  GPE(I )+AR*VEL 
WRITE(6,9000)  GRE (  I ) , G  I  E( I ) , TEMP 
IF ( I .EC. ITEMP.AND.MARK.EQ.2)  WRITE (6 f 900 l) 
55  CONTINUE 
C 
C  PUT  EIGENVALUES  INTO  SINGLE  PRECISION  VECTORS  TO  PASS  TO 
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C    SUBROUTINE    TO    DO    PLOTTING    FOR    ODD    FUNCTIONS. 
C 

DO    6  0    1  =  1  tN 

GR4(I)    =    SNGLCGRO(I)) 
6C    GI4( I)     =    SNGL(GIOd) ) 

WPITE(6,90G4) 

CALL    FLCTP(GR4,GI4,N,0) 

WRITE(6,9005)     N,  RE Y, A, B, VEL 
C 

C    SIMILARLY    PLOT    EIGENVALUES    FOR    EVEN    EI GENFUNCT IONS. 
C 

DO   65    1  =  1, N 

GR4(  I)    =    SNGL(GREU)  ) 
65    GI4( I)    =    SNGL(GIE( I) ) 

HRITE(6,9004) 

CALL    PLCTP<GR4,GI4,N,0J 

WRITE(6,9005)     N, RE Y, A, B, VEL 
C 

GO    TO    1 
IOC    WRITE(6,9004) 

STOP 

END 
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c 

C  PRCGRAM  #1A 

C 

C  PROGRAM  TO  PRINT  EIGENVALUES 

C  FOR  THE  3-D  POISEUILLE  FLOW  PROBLEM 

C  AND  PLOT  CN  CALGCMP  PLOTTER 

C 

c 

c 

c 

C      THIS  PROGRAM  SOLVES  THE  LINEARIZED  NAVIER  STOKES 

C      EGUATICN  FOR  POI  SEUiuLE  FLOrJ.   THE  EIGENVALUES 

C      RESULTING  FROM  THE  FINITE  DIFFERENCE  APPROXIMATION 

C      ARE  PRINTED. 

C 

C      INPUT 

C 

C  THE  FOLLOWING  MUST  BE  INPUT  TO  THE  PROGRAM  AS  DATA 

C  USING  NAMELIST,  'LIST*  . 

C 

C  N  -  HALF  OF    THE  NUMBER  OF  FINITE  DIFFERENCE  GRID 

C  POINTS  ACROSS  THE  CHANNEL  NOT  INCLUDING  THE  END 

C  POINTS.   N  MUST  BE  .LE.  MDIM,  WHICH  IS  THE 

C  DIMENSION  OF  THE  MATRICES  IN  THIS  PROGRAM. 

C  DEFAULTED  TO  THE  VALUE  OF  MDIM,  THAT  IS,  THE 

C  DIMENSION  OF  THE  MATRICES.   SEE  PROGRAM  BELOW  FOR 

C  THE  DEFAULT  VALUES. 

C 

C  REY  -  THE  REYNOLDS  NUMBER   <REAL*8) 

C  DEFAULT  VALUE  =  6000. 

C 

C  AR,AI,BR,BI  -  THE  REAL  AND  IMAGINARY  PARTS  OF  THE 

C  PERTURBATION  WAVE  NUMBERS   (REALMS).   DEFAULTED 

C  TO  0.,  1.,  0.  AND  0.  RESPECTIVELY. 

C 

C  VEL  -  THE  VELOCITY  OF  THE  MOVING  COORDINATE 

C  REFERENCE  SYSTEM  FOR  WHICH  THE  STABILITY  IS 

C  DETERMINED.   (REAL*8).   DEFAULTED  TO  0. 

C 

C  GRSCL,GISCL  -  THE  X  AND  Y  SCALE  OF  THE  OUTPUT 

C  GRAPH  IN  UNITS  PER  INCH.   BOTH  DEFAULTED  TO  0.2. 

C  SEE  SUBROUTINE  DRAW  FOR  RESTRICTIONS  ON  THEIR 

C  VALUES 

C 

C  NW,NH    -    THE    WIDTH    AND    HEIGHT    OF    THE    OUTPUT    GRAPH 

C  IN    INCHES.       DEFAULTED    TO    6    AND    8    RESPECTIVELY. 

C  SEE    SUBROUTINE    DRAW    FOR    RESTRICTIONS    ON    THEIR 

C  VALUES 

C 

C  TITLE  -  A  VECTOR  OF  LENGTH  12  WHICH  IS  PRINTED 

C  WITH  THE  GRAPH  (REAL*8) 

C 

C      OUTPUT 

c 

C  THE  OUTPUT  OF  THIS  PROGRAM  IS  A  TABULATION  OF  THE 

C  EIGENVALUES.   TWO  LISTS  ARE  PRINTED,  ONE  FOR  THE 

C  EIGENVALUES  CORRESPONDING  TO  EVEN  EI GENFUNCT IONS 

C  ANC  GNE  FOR  THOSE  CORRESPONDING  TO  ODD  EIGEN- 

C  FUNCTIONS.   THE  STABILITY  OF  EACH  EIGENVALUE  IS 

C  PRINTED  AND  THE  LEAST  STABLE  EIGENVALUE  IS  MARKED 

C  WITH  ASTERISKS.   THE  EIGENVALUES  ARE  GRAPHED  BY 

C  THE  CALCCMP  PLOTTER. 

C 

C      SUBROUTINES 

C 

C  THIS  PROGRAM  CALLS  THE  SUBROUTINE  'DEIGEO1  TO 

C  SOLVE  FOR  THE  EIGENVALUES  AND  SUBROUTINE  'DRAW' 

C  TO  GRAPH  THE  RESULTS. 

C 

c 

c 

c 
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IMPLICIT  REAL*8  (A-H,0-Z) 

CCMPLEX*16  A,B 

REAL**  GR4(30) ,GI4(30) 

REAL* 8  GRE(60) ,GIE(60) ,GR0(6  0) ,GI0(60) 

C0MPLEX*16  XMAT(60,30,3) 

CCMPLEX*16  YMAT(60,30),WVEC(60) ,BMAT( 5,60) 

EQUIVALENCE (YMAT(1,  1)  ,XMAT(  1,1,3)), 

*  (BMAT( 1,1) ,XMAT( 1,1,3)), 

*  (WVECt 1)  ,XMATi 1,6,3)  ) 
REAL*8  LABEL, TITLE(12) 

REAL*4  GRSCL,GISCL 

DATA  LABEL/1  »/ 

DATA  TITLE/'HARRISCN* , •   0279   »,10*»  •/ 

NAMELIST  /  LIST  /  N, REY, AR , AI , BR ,b I , VEL, GRSCL , GI SCL, 

*  NW,NH, TITLE 
C 

C  INITIALIZE  VARIABLES  (SET  DEFAULT  VALUES) 
C 

MDIM  =  60 

N  =  60 

REY  =  6C00D0 

AR  =  0CO 

AI  =  ICO 

BR  =  OCO 

BI  =  OCO 

VEL  =  ODO 

GRSCL  =  .20 

GISCL  =  .20 

NW  =  6 

f\H  =  8 
C 

C  READ  NAMELIST  AND  SET  ALPHA  AND  BETA 
C 

1  READ(5,LIST,END=100) 

A  =  DCNFLXUR,  AI  ) 

B  =  DCfPLX(BR,BI  ) 
C 

C  PRINT  INPUT  VALUES  AS  PAGE  HEADING  FOR  EIGENVALUE  LIST. 
C 

KRITE( 6,9004) 

9004  FCRMAT  CI') 

V»RITE(  6,9005)  N,  RE  Y,  A,  B,  VE  L 

9005  FCRMATf  N  =,,I4,/,'  REY  =  •  ,  F  10.  2,  8X  ,  ■  ALPH  A  =  •  , 

*  2F12.7,8X, 'BETA  = • , 2F1 2 . 7 , /  ,  •  VEL  =  ',F7.2) 
C 

C  CALL  SLBFCLTINE  TO  SOLVE  FOR  EIGENVALUES 
C 

CALL  DEIGEO(A,B,REY,N,MDIM,GRE,GIE,GRO,GIO,XMAT, YMAT, 

*  BMAT,WVEC) 
C 

C  DETERMINE  LEAST  STABLE  EIGENVALUE 
C 

TEMP  =  -1D10 

MARK  =  1 
C 

CO  20  1  =  1, N 

IF(GFO(I )+AR*VEL.LT.TEMP)  GO  TO  20 

TEMP  =  GRO(I )+AR*VEL 

ITEMP  =  I 
20  CONTINUE 
C 

CO  40  1  =  1, N 

IF(GREU)+Ak*VEL.LT.TEMP)  GO  TO  40 

TEMP  =  GRE(I )+AR*VEL 

ITEMP  =  I 

MARK  =  2 
40  CONTINUE 
C 
C  LIST  EIGENVALUES  FOR  ODD  EIGENVECTORS 


C 


WRITE(6,9003) 
9003  FCRMAT(///,6X, 'GAMMA  REAL1  , 5X, ■ GAMMA  I  MAG*  ,12X , ■ ST AB« ) 
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WRITE(6,9006) 
90C6  FGRMATC • OE IGENVALUES  FOR  ODD  EIGENVECTORS',/) 
C 

CC  50  1  =  1,  N 
TEMP  =  GPOU  )+AR*VEL 
WRITE<  6,9  000)  GRC(  I ) ,G  10 ( I ) ,TEMP 
IF  I  I  .EG.  I  TEMP. AND. MARK. EQ.l)  WRITE (6,9001) 
50  CCNTINLE 
90C0  FORMAT (  »0 •  ,  1P2D1 5. 4, 1PD2  0. 4  J 
9001  FORMAT (• +  *  ,52X, •***« ) 
C 

C  LIST  EIGENVALUES  FOR  EVEN  EIGENVECTORS. 
C 

WPITE(6,9007) 
90C7  FOPMATCOEIGENVALUES  FOR  EVEN  EIGENVECTORS*,/) 
DC  55  1=1, N 
TEMP  =  GRE(I)+AR*VEL 
WRITE(6,9000)  GRE {  I )  ,G  I  E ( I ) , TEMP 
IF (I .EC. I  TEMP. AND. MARK. EQ. 2)  WRI TE (6 ,900 1 ) 
55  CCNTINLE 
C 

C  PUT  EIGENVALUES  INTO  REAL*4  VECTOR  IN  GROUPS  OF  30  OR  LESS 
C  TO  PASS  TO  SUBROUTINE  DRAW.   CALL  DRAW  REPEATEDLY  UNTIL 
C  ALL  EIGENVALUES  FOR  ODD  EIGENVECTORS  ARE  PLOTTED.   DO  THE 
C  SAME  FOR  THE  EIGENVALUES  FOR  EVEN  EIGENVECTORS. 


C 


NGCOCE  =  1 
IE  =  30 
IMODE  =  0 
IF130.GT.N)  IE  =  N 

60  DO  65  1=1, IE 

GR4(I)  =  SNGL(GRC(IM0DE*30*IJ) 
65  GI4(I)  =  SNGLCGICt IM0DE*30+I)) 

CALL  CRAW( IE,GR4,GI4,NGC0DE,3,LABEL,TITLE, 

*  GRSCL,GISCL,NH,NW,2,2,NW,NH,0, LAST) 

NGCODE  =  2 

IF((IMCDE*1)*30+1.GT.N)  NGCODE  =  3 
DO  7  0  1  =  1,  IE 

GR4(I)  =  SNGL(GRE(  1M0DE*30+I )) 
7C  GI4(I)  =  SNGL(GIE(IMCDE*30  +  D) 

CALL  CRAW( IE,GR4,GI4,NGC00E,1,LABEL,TITLE, 

*  GRSCL,GISCL,NH,NW,2,2,NW,Nh,0,LAST) 

IF(NGCCCE.EQ.3)  GO  TO  80 

IMCDE  =  IMCDE  +  1 

IF(( IMCDE+li*30.GT.N)  IE  =  N  -  IM0DE*30 

GO  TO  60 

80  GO  TO  1 
100  STCP 
END 


131 


c 

C  PRCGRAM  #2 

C 

C  PROGRAM  TO  FIND  THE  SHAPE  OF 

C  THE  STABILITY  CURVES 

C  FOR  3-D  POISEUILLE  FLOW 

C 

c 

c 

C      THIS  PROGRAM  FINDS  THE  SHAPE  OF  THE  STABILITY  CURVES. 

C      IT  PERFORMS  FIVE  FUNCTIONS  DEPENDING  ON  THE  VALUE  OF 

C      THE  INPUT  PARAMETER  'MODE1.   INPUT  IS  BY  THE 

C      NAMELIST  'LIST'.   OTHER  INPUT  AND  OUTPUT  ARE 

C      DETERMINED  BY  THE  VALUE  OF  MODE  AS  WELL  AS  THE 

C      FUNCTION  PERFORMED.   IN  ALL  CASES,  ONLY  THE  PARAMETERS 

C      AI  AND  REY  ARE  VARIED.   OTHER  PARAMETERS  ARE  ALWAYS 

C      hELD  CONSTANT  THROUGHOUT  THE  CALCULATIONS. 

C 

C      MODE  =  1  

C 

C      THE  STABILITY  AT  ANY  GIVEN  POINT  IS  FOUND 

C 

C      INPUT 

C 

C  AR,AI,BR,BI  -  THE  REAL  AND  IMAGINARY  PARTS 

C  OF  ALPHA  AND  BETA.  <REAL*8) 

C 

C  VEL  -  THE  VELOCITY  OF  THE  MOVING  COORDINATE 

C  FRAME  FOR  WHICH  THE  STABILITY  IS  TO  BE 

C  DETERMINED.  (REAL-8) 

C 

C  REY  -  THE  REYNOLDS  NUMBER.  (REAL*8) 

C 

C      OLTPUT 

C 

C  THE  PROGRAM  PRINTS  THE  INPUT  AND  THE  RESULTING 

C  STABILITY. 

C 

C      MODE  =  2 

C 

C      FOR  A  GIVEN  VALUE  OF  REY,  AI  IS  ADJUSTED  TO  FIND  THE 

C      POINT  WHICH  HAS  THE  DESIRED  STABILITY. 

C 

C      INPUT 

C 

C  AR,BR,BI,REY,VEL  -  VALUES  TO  BE  HELD  CONSTANT 

C  DUPING  THE  SEARCH  FOR  THE  VALUE  OF  AI. 

C 

C  AI  -  AN  INITIAL  GUESS. 

C 

C  CAI  -  AN  ESTIMATE  OF  THE  ERROR  IN  THE  GUESS 

C  FOR  AI.  (REAL*8) 

C 

C  STAB  -  THE  STABILITY  FOR  WHICH  THE  VALUE  OF  AI  IS 

C  CESIREC.  (REALMS) 

C 

C      OUTPUT 

C 

C  THE  INPUT  VALUES  INCLUDING  THE  INITIAL  GUESS  FOR 

C  AI  ARE  PRINTED  ALONG  WITH  THE  FINAL  APPROXIMATION 

C  FOR  AI. 

C 

C      MODE  =  3  

C 

C      FOR  A  GIVEN  VALUE  OF  AI,  REY  IS  ADJUSTED  TO  FIND  THE 

C      POINT  WHICH  HAS  THE  DESIRED  STABILITY. 

C 

C      INPUT 

C 

C  AR,AI,BR,BI  -  VALUES  TO  BE  HELD  CO^TANT  DURING 

C  ThE  SEARCH  FOR  THE  VALUE  OF  REY. 

C 
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C  REY  -  AN  INITIAL  GUESS. 

C 

C  DREY  -  AN  ESTIMATE  OF  THE  ERROR  IN  THE  GUESS  FOR 

C  REY.  (REAL*8) 

C 

C  STAB  -  THE  STABILITY  FOR  WHICH  THE  VALUE  OF  REY 

C  IS  CESIREO. 

C 

C      OUTPUT 

C 

C  ThE  INPUT  VALUES  INCLUDING  THE  INITIAL  GUESS 

C  FOR  REY  ARE  PRINTED  ALONG  WITH  THE  FINAL 

C  APPROXIMATION  FOR  REY. 

C 

C      MODE  =  4  

C 

C      THE  POINT  CF  MINIMUM  REY  ON  THE  NEUTRAL  STABILITY 

C      CURVE  IS  DETERMINED. 

C 

C      INPUT 

C 

C  AP,ER,BI,VEL  -  VALUES  TO  BE  HELD  CONSTANT  DURING 

C  THE  DETERMINATION  OF  THE  MINIMUM  REY. 

C 

C  AI,PEY  -  INITIAL  GUESS  FOR  THE  POINT  OF 

C  MINIMUM  REY. 

C 

C  DAI, DREY  -  AN  ESTIMATE  OF  THE  ERROR  IN  THE 

C  GUESSES  FOR  AI  AND  REY. 

C 

C      OUTPUT 

C 

C  ThE  PROGRAM  PRINTS  THE  INITIAL  GUESS,  THE 

C  SUCCESSIVE  APPROXIMATIONS  AND  THE  FINAL 

C  ESTIMATE  FOR  THE  VALUES  OF  AI  AND  REY. 

C 

C      MODE  =  5  

C 

C      TF.E  POINT  CF  MAXIMUM  INSTABILITY  IS  DETERMINED. 

C 

C      INPUT 

C 

C  AR,BR,BI,VEL    -    VALUES    TO    BE    HELD    CONSTANT 

C 

C  AI,REY    -     INITIAL    GUESS     FOR    THE    POINT    OF    MAXIMUM 

C  INSTABILITY. 

C 

C  CAIfDREY  -  AN  ESTIMATE  OF  THE  ERROR  IN  THE 

C  GUESSES  FOR  AI  AND  REY. 

C 

C      OUTPUT 

C 

C  THE  PROGRAM  PRINTS  THE  INITIAL  GUESS,  THE 

C  SUCCESSIVE  APPROXIMATIONS  AND  THE  FINAL  ESTIMATE 

C  FOR  ThE  VALUES  OF  AI  AND  REY. 

C 

C       

c 

C      SUBROUTINES  REQUIRED 

C 

C  THIS  PROGRAM  CALLS  THE  SUBROUTINE  ,DFIT2f  AND  USES 

C  THE  FUNCTION  'EIGLS'. 

C 

c 


c 


IMPLICIT  REAL*8  (A-H.O-Z) 

NAMELIST  /  LIST  /  REY, AR ,AI , BI , DAI , DREY, VEL, BR, MODE , 


STAB 
WRITE (6,9000) 
900C  FORMAT (« 1«) 
C 
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C  INITIALIZE  VARIABLES 
C 

STAB  =  CDO 

VEL=ODO 

BR  =  OCO 

BI  =  OCO 

AR  =  OCO 

AI  =  ICO 

CAI  =  .OlOCO 

DREY  =  250D0 

MCDE  =  1 
C 

C  READ  NAMELIST  VARIABLES 
C 

1  READ(5,LIST,END=100) 

PEY1  =  REY 

AITEMP  =  AI 
C 

C  JUfP  TC  TYPE  OF  SOLUTION  CALLED  FOR  BY  VALUE  OF  MODE 
C 

GC  TO  (10,20,30,40,50)  ,MODE 
C 

C      MODE  =  I 

C  FIND    STAB    FOR    SINGLE    POINT 

C 

10    WRITE     (6,9100)     REY , AR, AI , BR , BI , VEL 

9100  FORMAT (///,21X,« STABILITY    AT    POI NT . . . ■ ,/  ,2 IX , «REY    =  «, 

*  F8.lt/t2.lXt'  AR    =  •  ,F12.5,5X, ' AI    =  «  ,  Fl  2.  5,  /  ,  21X  , 

*  »BR    =• ,F12.5,5X, »BI    =  « ,F 12. 5 ,/ , 21X , • FOR    VEL    =  ', 

*  F6.2) 
C 

C    USE    EIGLS    TO    FIND    STABILITY 
C 

Gl    =    EIGLS(REY,AR,AI ,BR,BI ,VEL) 

WRITE(6,9I01)    Gl 

9101  FCPMAT (/,21X,« STAB  =«,1PD12.4) 
GC  TO  1 

C 

C      MODE  =  2..... 

C  FIND    AI    GIVEN    REY    AND    STAB 

C 

2C    WRITE(6,9200)     STAB  ,REY  ,  AI , BR ,B I , AR , VEL ,D AI 
920C    FGPMAT(///,21X,'F0R    STAB    = • , F9 . 5 ,5X, ■ AT    REY    =,,F8.1, 

*  /,21X, 'INITIAL    GUESS.. .       AI    =•  ,F 12 . 5 ,/  ,21X , 

*  «GIVEN    BR    =• ,F12.5,5X, 'BI    = • ,F 12 . 5 ,/ ,2 7X , • AR    =«, 

*  F12.5,5X,'VEL    =•  ,F6. 2, / , 21X , «D AI    =',F10.5,/) 
C 

C    CALCULATE    STABILITY    FOR     INITIAL    GUESS. 
C 

Gl    =    EIGLS(REY,AR,AI  ,BR,BI ,VEL) 

WRITE(6,9201)     AI,G1 

9201  FGFMAT(21X,'F0R    AI    = • , F9 .5 , 5X, • STAB    =»,1PD12.4) 
C 

C    ADD    CAI     TO    INITIAL    AI     AND    DETERMINE    STABILITY. 
C 

All    =    AI    +    DAI 

G2    =    EIGLS(REY,AR, AI 1,BR,BI ,VEL) 

WRITE(6,S201)    AI1,G2 
C 

C    USE    LINEAR    APPROXI MATI CN    TO    FIND    NEW    VALUE    OF    AI . 
C    CALCULATE    STABILITY. 
C 

AI2    =    All    +    <AI1-AI)*(STAB-G2)/(G2-G1) 

G3    =    EIGLS(REY,AR,AI 2,BR,BI ,VEL) 

WRITE(6,9201)     AI2,G3 
C 

C  USE  SUBROUTINE  DFIT2  TO  OBTAIN  SECCND  ORDER  APPROXIMATION 
C  FOR  FINAL  VALUE  OF  AI. 
C 

CALL  DFIT2(AI , A I  1, AI 2 , Gl ,G2 , G3 , STAB, AI  f :, DUMMY,DUMY2) 

hRITE(6f52C2)  AIF 

9202  F0RMAT(/,21X, 'FINAL  ESTIMATE...  AI  =«,F12.5) 
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GC    TO    1 

C  MODE    =    3 

C  FIND    REY    GIVEN    AI    AND    STAB 

C 

30    WRITE(6,9300)     STAB , A  I , RE Y, BR, BI , AR , VEL ,DREY 

9300  FCRMAT(///t21X,' FOR    STAB    = ' , F9. 5 , 5X, • AT    AI     =«,F12.5, 

*  /,21X, ' INITIAL    GUESS...       REY    =  •  ,  F8 . 1 , / , 21X , 

*  'GIVEN    BR    -• fF12.5,5X, «BI    =  '  ,F 12 . 5 ,/ ,27X , 

*  «AR    =  •  ,F12.5,5X, 'VEL    =  '  ,  F6  .2 ,/ , 2 IX , • DREY    =•, 

*  F3.1,/) 
C 

C    CALCULATE    STABILITY    FOR    INITIAL    GUESS 
C 

Gl    =    EIGLS(REY,AR,AI ,BR,BI ,VEL) 

URITE(6,9301 )  REY,G1 

9301  FCPMAT(21X,'F0R    REY    =  '  ,F8. 1 , 5X i • ST AB    =',1PD12.4) 
C 

C    ACC    CREY    TC     INITIAL    REY    AND    DETERMINE    STABILITY. 
C 

REY1    =    REY    +    DREY 

G2    =    EIGLS<REY1,AR,AI,BR,BI,VEL) 

WF1TE(6,9301)     REY1.G2 
C 

C    USE    LINEAR    APPROXIMATION    TO    FIND    NEW   VALUE    OF    REY. 
C    CALCULATE     STABILITY. 
C 

REY2    =    REY1    +    ( REY 1-REY ) *( STAB-G2 ) / ( G2-G 1 ) 

G3  =  EIGLS(REY2, AR , AI , BR , B I , VEL) 

WRITE(6,9301)  REY2,G3 
C 

C  USE  SUBROUTINE  DFIT2  TO  OBTAIN  SECCND  ORDER  APPROXIMATION 
C  FOR  FINAL  VALUE  OF  REY. 
C 

CALL  DFIT2(REY,REY1,REY2,G1,G2,G3, ST AB , REY F , DUMMY , 
.  *     DUMY2) 

fcRITE(6,9302)  REYF 

9302  F0FMAT(/,21X, 'FINAL  ESTIMATE...   REY  =',F8.1) 
GG  TO  1 

C 

C MODE  =  4 

C         FIND  POINT  OF  MINIMUM  REYNOLDS  NUMBER 

C  ON  THE  NEUTRAL  STABILITY  CURVE 

C 

C  OUTPUT  VALUES  FOR  INITIAL  GUESS 

C 

40    WRITE(6,9400)    AI  ,  REY  ,0 AI ,DREY , B I  , A R , VEL , BR 

9400  FCRMAT( • 1' ,////, 21X, 'POINT    UF     MINIMUM    REY    ON    NEUTRAL', 

*  '     STABILITY    CURVE'  ,/ ,21X, 

*  'INITIAL    GUESS...       AI    =',F12.5,5X, 

*  'REY1    =' ,F8. 1,/,21X, 'DAI     =  ' ,  Fl  2.  5,  5X  , '  DREY    =' , 

*  F8.1,/,21X,'B1    =• ,F12.5,5X,' AR    =',F12.5, 

*  /,21X,'VEL    ='  ,F6.2,5X, 'BR    =  '  ,F  12  .5  ,/// ) 
C 

C    FOR     INITIAL    REYNOLDS    NUMBER,    CHOOSE    TW C    MORE    VALUES    FOR 
C    AI    ANC    USE    PARABOLIC    CURVE    FIT    TO    FIND    MINIMUM    STABILITY 
C    WITH    CORRESFCNCING    VALUE    FOR    AI 
C 

AIM    =    AI-CAI 

AIP    =    AH-DAI 

Gl    =    EIGLS(REY1, AR , A IM , BR, BI , V  EL ) 

WRITE(6,94C1)    AIM,G1 

9401  F0RMAT(21X,'REY1.. .  AI    = • , F  12. 5 , 5X , • FOR    STAB    =', 

*  1PC12.4) 

G2    =    EIGLS(REY1,AR,AI,BR,BI,VEL) 

WRITE(6,9401)     AI,G2 

G3  =  EIGLS(REY1, AR,AIP,BR,BI ,VEL) 

URITE(6,94C1 )  AIP,G3 

CALL  DFIT2(A!M,AI,AIP,G1 ,G2 , G3 , 1D2 , DUMMY ,AI1,GMIN1) 

WR I T  E ( fc , 9  4C2  >  All, GM I N 1 

9402  FORMATCO1  f20Xf 'AI  =  *  ,  Fl  2.  5  ,  5X  ,  '  FOR  MAX  G  ='  , 

*  1PC15.4,/) 
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C    CHCCSE    NEXT    GUESS    FOR    REYNOLDS    NUMBER    USING    DREY 
C 

REY2    =    REY1+DREY 

IF(GMIM  .GT.ODOJ    REY2    =    REY1-DREY 
AI    =    All 

V»RITE(6,9403)    AI.REY2 
9403    F0RMAT(///f21X,» SECOND    GUESS...       AI    =«,F12.5,5X, 
*  'REY2    =' ,F8.1,/) 


C 
C 
C 
C 


C 

c 
c 


c 
c 
c 
c 


c 
c 

c 


c 
c 
c 
c 
c 


FIND    MINIMUM    STABILITY    AND    CORRESPONDING    AI     FOR    SECOND 
GUESS    OF    REY 

AIM    =    AI-DAI 

AIF    =    AI+CAI 

Gl    =    EIGLS(REY2,AR ,  A IM ,BR, BI , VEL ) 

fcRITE(6,94G4)  AIM.G1 

9404  FCFMAT(21X,»  REY2..  .      AI  =  ■  ,  Fl  2.  5  ,  5X  ,  •  FOR  STAB  =  «, 

*  1PD12.4) 

G2    =    EIGLS(REY2,AR,AI,BR,BI ,VEL) 

WRITE(6,9404)    AI,G2 

G3    =    EIGLS(REY2,AR,AIP,BRtBIfVEL) 

WRITt(6,9404)     AIP,G3 

CALL    0FIT2(AIM,AI,AIP,G1,G2,G3 ,1D2, DUMMY ,AI2 ,GMIN2) 

WRITE(6,9402)    AI2,GMIN2 

USE    LINEAR    FIT    TO    DETERMINE    THIRD    GUESS    FOR    REY    AND    AI 

REY3    =     REYl-( REY 1-REY2 )*GMI Nl/ (GMI Nl-GMI  N2 ) 
AI    =    AI1-(AI1-AI2)*GMIN1/<GMIN1-GMIN2> 
VnRITE(6,9405)     AI,REY3 

9405  FORMAT!///, 21X,« THIRD    GUESS...       AI     =',F12.5,5X, 

*  *REY3    =  » ,F8.1,/) 

FEFEAT    PFOCECURE    TO    FIND    MINIMUM    STABILITY    AT    THIRD    GUESS 
FCP    PEY 

AIM    =    AI-CAI 

AIP    =    AI+CAI 

Gl    =     EIGLS(REY3, AR,AIM,BR,BI ,VEL) 

hRITE(6,9406)  AIM,G1 

9406  FGRMAT(21X, •REY3...      AI  = • , F 12. 5, 5X , • FOR  STAB  =  ', 

*  1PC12.4) 

G2    =    EIGLS(REY3, AR , AI , BR , BI , VEL ) 

hRITE(6,9406)     AI ,G2 

G3    =    EIGLS(REY3,AR,AIP,BR,BI ,VEL) 

hRITEC 6,9406)  AIP,G3 

CALL  DFIT2(AIM,AI , A I P , Gl ,G2 , G3  , 1 D2 , DUMMY ,AI3,GMIN3) 

WRITE(6,94G2)  AI3,GMIN3 

USE    PARABOLIC    FIT    TO    DETERMINE    FINAL    GUESS    FOR    REY    AND    AI 

CALL    DFIT2(AI1,AI2,AI3,GMIN1,GMIN2,GMIN3 ,ODO,AIF, 

*  DUMMY,DUMY2J 

CALL    DFIT2(REY1,REY2,REY3,GMIN1,GMIN2,GMIN3,0D0,REYF, 

*  DUMMY,DUMY2) 
WRITE(6,94C7)  AIF,REYF 

94C7    FORMAT!///, 21X,» FINAL    ESTIMATE...       AI    =»,F12.5,5X, 

*  »REY    =*,F8.1, ////////) 
AI    =    AITEMF 
GO    TO     1 


MODE  =  5 

FIND  POINT  OF  MINIMUM  STABILITY 
(THAT  IS,  FOR  MAX  VALUE  OF  STAB) 

50  WRITE(6,9500)  A  I  ,  REY  ,DAI ,DREY , BI , A R, VEL, BR 
950C  FORMAT!  • 1«  ,///, 20X  ,« POINT  OF  MINIMUM  STABI LI TY« , / , 2 IX , 

*  'INITIAL  GUESS...   A  I  = • , F12  .5 , 5X,  'REY  1  =',F3.1, 

*  /,21X,*DAI  =» ,F12.5,5X,,DREY  =,,F8.1,/,21X, 

*  'BI  =•  ,P12.5,5X,  fAR  =• ,F12.5,/,21X,' VEL  =f,F6.2, 

*  5X, •BR  =*  ,F12.5,///) 
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AIM    =    AI    -    DAI 

AIP    =    AI    ♦    DAI 

Gl    =     EICLS(REY1,AR  ,AIM,BR,BI ,VEL) 

URITE(6,95C1)     AIM.G1 


9501  FCFMAT(21X,'REY1 
*     1PD12.4) 
G2  =  EIGLSCREYli 
hRITE(6,9501J  AI 

G3  =  EIGLS(REY1, 
hRITE(6,95Cl )  AI 


AI  =• ,F12.5t5X, 'FOR  STAB  =  », 


AR,AI ,BR,BI ,VEL) 
G  ? 

arTaip,br,bi,ved 

P,G3 
CALL  DFIT2CAIMf AIfAIPfGl»G2tG3tlD2 

WRITE(6,9  5  02)    AI1,GMIN1 


,DUMMY,AI1 ,GMIN1) 
FOR  *0»,20X,'AI    =«,F12.5  :0P    STAB    = • , 1 PD15 .4 , / ) 


95 


REY2  =  PEY1  +  CREY 
AI  =  All 

kRITE(6,9503)  AI ,REY2 
03  FORMATI///^^,1  SECOND 
*     'REY2  =  • fF8.li/) 


GUESS 


AI  =«,F12.5,5X 


AIM  =  AI  -  DAI 
AIP  =  AI  +  DAI 
Gl  =  EIGLS(REY2f ARfAIMtBRfBIfVEL) 

V\RITE(6,9504)  AIM.G1 
9504  FCPMAT(21X,,REY2...     AI  = • , Fl 2 .5 ,5X , ■ FOR  STAB  =  «, 
*     1PD12.41 

G2    =    EIGLS(REY2,AR,AI,BR,BI ,VEL) 

WRITE(6,9504)    AI,G2 

G3    =    EIGLS<REY2, AR,AIP,BR,BI ,VEL) 

WRITE(6,9504)     AIF,G3 

CALL  UF I T2( AIM, A I , A I P , Gl ,G2 ,G3 , 1D2 , DUMMY ,AI2,GMIN2) 

WRITE(6,9502)  AI2,GMIN2 


95C5 


REY3    =    REY1    -    DREY 

IF(GMIN2.GT.GMIN1)    REY3    =    REY2    +    D 
AI    =    AI2    *    (All-AI2)*(REY3-REY2)/l 
WRITE(6,9505 J    AI,REY3 
FaPMATl///,21X,« THIRD    GUESS...       AI 
'REY3    =' ,F8.1,/) 


REY 
REY1-REY2) 

=  •  ,F1 2.5, 5X, 


9506 


C 
C 


AIM    =    AI    -    DAI 

AIP    =    AI    +    DAI 

Gl    =    ElGLS(REY3,ARfAIMfBRf BlfVEL) 

WRITE (6,9  506)     AIM,G1 

FCPMAT(21X,« REY3...  AI    =  ■ , Fl 2 .5 ,5X, • FOR    STAB    =  ', 

1PD12.4J 
G2    =    EIGLS(REY3,AR,AI,BR,3I ,VEL) 
hRITE(6,9506)    AI,G2 
G3    =     EIGLS(REY3,AR,AIP,BR,BI  ,VEL) 
WRITE(6,9506)     AIP,G3 

CALL    DFIT2(AIMiAI ,AIP,G1 ,G2 , G3 , 1D2 , DUMMY ,AI3,GMIN3J 
WRITE(6,9502)    AI3,GMIN3 

CALL    DFIT2(REY1,REY2,REY3,GMIN1,GMIN2,GMIN3,1D0,DUMMY, 
REYF,GF) 


A    =     ( ( AIl-AI2)/(REYl-REY2)-( AI1-AI 

*  /( (REY1#*2-REY2**2)/(REY1-REY2 

*  -(REYl**2-REY3**2)/< RE Yl- REY 3) 
B  =  (AI1-AI2-A=MREY1**2-REY2**2)  )/ 
C    =    AI1-A*REY1**2-B*REY1 

AIF    =    A*PEYF**2+B*REYF+C 


3)/(REYl-REY3)) 

J 

) 

(REY1-REY2) 


WRITE(6,9507)     AIF,REYF,GF 
95C7    FGRMAT(///,21X, 'FINAL    ESTIMATE...       AI     =•  ,F12.5,5X, 
*  'REY    =«fF8.1t/i41Xf «STAB    =', 1PD1 5.4, //////// ) 

GO    TO    1 

IOC    VsRITE(6,9000) 
STOP 
END 
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C FUNCTION  EIGLS , 

C 

C      PUPPCSE 

c 

C        EIGLS  RETURNS  THE  STABILITY  OF  THE  LEAST  STABLE 

C        EIGENVALUE  FUR  PLANE  POISEULLE  FLOW  AS  DETERMINED 

C        FBCM  A  REFERENCE  FRAME  MOVING  IN  THE  X  DIRECTION 

C        WITH  VELCCITY  VEL. 

C 

C      USAGE 

C 

C        STAB  =  EIGLS(REY,AR,AI,BR,BI,VEL) 

C 

C      CESCRIPTICN  OF  PARAMETERS 

C 

C        THE  FOLLOWING  MUST  BE  SET  BY  THE  CALLING  PROGRAM.. 

C  REY,AR,AI,BR,BI,VEL 

C 

C        REY  -  THE  REYNOLDS  NUMBER  (REAL*3) 

C        AR,AI,BR,BI  -  THE  REAL  AND  IMAGINARY  PARTS  OF  THE 

C  PERTURBATION  WAVE  NUMBERS,  ALPHA  AND  BETA 

C  (REAL*8J 

C 

C        VEL  -  THE  VELOCITY  OF  THE  COORDINATE  SYSTEM  FOR 

C  WHICH  THE  STABILITY  IS  DESIRED.   (THE  MINIMUM, 

C  AVERAGE,  AND  MAXIMUM  VELOCITIES  OF  THE  FLUID  ARE 

C  O.O,  1.0,  AND  1.5  RESPECTIVELY)   (REAL*8) 

C 

C  THE    VALUE    OF    THE    STABILITY    OF    THE    LEAST    STABLE 

C  EIGENVALUE     IS    RETURNED    IN    'EIGLS'.       EIGLS    MUST    BE 

C  REAL*8    IN    THE    CALLING    PROGRAM. 

C 

C      OTHER  ROUTINES  NEEDED 

C 

C        DEIGEO 

C 

C 

C 

FUNCTION    EIGLS (REY,AR,AI ,BR,BI ,VEL) 

IMPLICIT    REALMS     (A-H,0-Z) 

REAL* 3    GRE(60) ,GIE(60J , GP0(60) ,GIO(60) 

CCMPLEX*16    XMAT(bO,30,3) 

COMPLEX* 16    YMAT(60,30) ,WVEC(60} ,BMAT( 5,60} 
E GUI  VALENCE (Y MAT (  1, 1)  ,XMAT( 1,1,3)  ) , 
*  (BMAT(lfl) ,XMAT< 1,1, 3) ), 


-»■  IDI'IHI    U  JlMA  I'lH   1  \    X  ,   X  ,    D 

*  (WVEC(l)  ,XMAT(1,6,3)  ) 


PDIM    =    60 

N    =    60 
C 

C  CALL  SUBROUTINE  TO  DETERMINE  EIGENVALUES  FOR  BOTH  EVEN 
C  AND  COC  EIGENF UNCTIONS. 
C 

CALL  CEIGEC(DCMPLX(AR, A I ), DCMPLX(BR,BI ) , REY , N, MD I M, 
*     GRE,GIE,GRC,GIO, XMAT , YMAT, BMAT ,WVEC) 
C 

C  DETERMINE  STABILITY  FOR  EACH  EIGENVALUE  SAVING  THE  VALUE 
C  FOP  THE  LEAST  STABLE  ONE. 
C 

TEMP  =  -1C10 

CO  10  1=1, N 

IF(GRE(  I  )+AR*VEL.LT.TEMP)  GO  TO  5 

TENP  =  GRE(I )+AR*VEL 
5  IF(GRG(I)+AR*VEL.LT.TEMP)  GG  TO  10 

TEMP  =  GRO(I )+AR*VEL 
10  CONTINUE 
C 

C  RETURN  THE  VALUE  OF  THE  MINIMUM  STABILITY  IN  •EIGLS1. 
C 

EIGLS  =  TEMP 

RETURN 

END 
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C FUNCTION  CHM1E1  AND  CHM2E1 

C  (CARTESIAN  COORDINATES) 
C 

c 

C  PURPOSE 

c 

C  CHM1E1  AND  CHiM2El  RETURN  THE  VALUES  (C0MPLEX*16)  OF 

C  ThE  COEFFICIENTS  FOR  THE  MATRICES  IN  THE  FINITE 

C  DIFFERENCE  FORM  OF  THE  LINEARIZED  NAVI tR-STOKES 

C  EQUATION  FOR  POISEUILLE  FLOW.   BOTH  FUNCTIONS  RESULT 

C  FRCM  THE  LINEAR  COMBINATION  OF  EQUATION  1  AND 

C  EQUATICN  3  TO  ELIMINATE  THE  VELJCITY  VECTOR 

C  POTENTIAL  COMPONENT  G  AND  ARBITRARILY  SETTING  THE 

C  CCMFCNENT  F  TO  ZERO.   SO,  THEY  ARE  THE  COEFFICIENTS 

C  FOR  THE  VECTOR  POTENTIAL  COMPONENT  H.   CHM2E1 

C  RETURNS  THE  TERMS  WHICH  ARE  COEFFICIENTS  OF  THE 

C  EIGENVALUE,  GAMMA,  AND  CHM1E1  RETURNS  THE  REMAINING 

C  TERMS. 

C 

C  LSAGE 

C 

C  XI  =  CHM1E1(K,Y) 

C  X2=CHM2E1(K,Y) 

C 

C  (CHM1E1  AND  CHM2E1  MUST  BE  DECLARED  C0MPLEX*16  IN 

C  CALLING  PROGRAM) 

C 

C  DESCRIPTION  OF  PARAMETERS 

C 

C  THE  FOLLOWING  PARAMETERS  MUST  BE  SET  BY  THE  CALLING 

C  PROGRAM 

C 

C  K  -  INDICATES  THE  POINT  ON  THE  FINITE  DIFFERENCE 

C  MESH.  RELATIVE  TO  ThE  CENTRAL  POINT  IN  THE  CENTRAL 

C  DIFFERENCING  SCHEME.   IF  THE  DIFFERENCE  IS  BEING 

C  FORMED  ABOUT  THE  N-TH  POINT  THEN  K=l  REFERS  TO 

C  THE  POINT  N-2,  K=2  REFERS  TO  THE  POINT  N-l, 

C  K=3  REFERS  TO  N,  K=4  REFERS  TO  N+l,  AND  K=5 

C  REFERS  TO  N+2. 

C  K  -  INDICATES  WHICH  POINT  ON  THE  FINITE  DIFFERENCE 

C  MESH  IS  REFERRED  TO  THE  CENTRAL  POINT.   IF  THE 

C  DIFFERENCE  IS  BEING  FORMED  ABOUT  THE  N-TH  POINT 

C  THEN  K=l  REFERS  TO  THE  POINT  N-2,  K=2  REFERS  TO 

C  THE  POINT  N-l,  K=3  REFERS  TO  N,  K=4  REFERS  TO  N+l, 

C  AND  K=5  REFERS  TO  N+2. 

C 

C  Y    -    THE    VALUE    OF    THE    POSITION    RELATIVE    TO    THE    CENTER 

C  OF    THE    CHANNEL.       THE    TWO    BOUNDARIES    ARE    AT    Y=+l 

C  AND    Y=-l. 

C 

C  CTFER  ROUTINES  NEEDED 

C 

C  NONE 

C 

c 

C 

c 


FLNCTICN    CHMlEKKfY) 

IMPLICIT    CCMPLEX*16    (A-H,C-Z) 

COMMON    /    COEFNT    /    A, B, G, REY, DEL 

REAL*8    REY, Y, DEL 
C 

C    THE    FOLLOWING    FUNCTIONS    (Ml)     EVALUATE    THE    COEFFICIENTS 
C    OF    THE    DERIVATIVES    OF    H    FOR    ALL    TERMS    EXCEPT    THOSE 
C    CONTAINING    GAMMA. 
C 

CH4MKY)     =    A/REY 

CH2MKY)    =    -1.5DO*A**2*(lDO-Y**2)+2DO*A*(A**2«-B**2) 

*  /REY 

CHOMHY)     =    -A*(( A**2+B**2)*( I.5DO*A*(lDu-Y**2)-<A**2 

*  +B**2)/REY)+3D0*A) 
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C  THE  REMAINING  FUNCTICNS  <M2)  EVALUATE  THE  COEFFICIENTS 
C  OF  THE  DERIVATIVES  OF  H  WHICH  ARE  ALSO  COEFFICIENTS 
C  CF  THE  EIGENVALUE  GAMMA. 
C 

CF2M2IY)  =  A 

CHOM21Y)  =  A*( A**2+B**2) 
C 

C  SET  LP  THE  FINITE  DIFFERENCE  VALUES  FOR  INDEX  K  FOR  Ml 
C 

GC  TC  (1,2,3,2,1 ),K 

1  CHM1E1    =    CH4M1(Y)/DEL**4 
GC    TC    1GO 

2  CHM1E1     =    -4D0*CH4M1( Y) /DEL**4+CH2M 1( Y ) /D EL**2 
GC    TO     100 

3  CHM1E1     =    6D0*CH4M1 (Y  )/ DEL**4-2D0*CH2M1 ( Y )/DEL**2 
*  +CHOMKY) 

ICC    RETURN 
C 
C    SET    LP    THE    FINITE    DIFFERENCE    VALUES    FOR    INDEX    K    FOR    M2 


C 


ENTRY    CHP2E1(K,Y) 

GC    TC    (11,  12,13,12,11)  ,K 

11  CHM2E1  =  (CDO,ODO) 
GC  TO  200 

12  CHM2E1  =  CH2M2(Y)/DEL**2 
GO  TO  200 

13  CHN2E1  =  -2D0*CH2M2(Y)/DEL**2+CH0M2(Y) 
2CC  FETUFN 

END 


140 


C SUBROUTINE  MULDBM 

C 

c 

C      FUPFCSE 

C 

C        MULDBN  PERFORMS  THE  MATRIX  MULTIPLICATION  BETWEEN  A 

C        SQUAPE  MATRIX  X  AND  A  BANGED  MATRIX  XB  WHICH  HAS 

C        BEEN  SET  UP  BY  SUBROUTINE  BMSET.   THE  RESULT  IS 

C        PLACEC  BACK  IN  X.   THAT  IS... 

C  X  =  (X) (XB) 

C 

C      LSAGE 

C 

C        CALL  NULDBM(X,XB,N,NB,MDIM,TEMPV) 

C 

C      CESCRIFTICN  OF  PARAMETERS 

C 

C  THE    FCLLOWING    MUST    BE    SET    BY    THE    CALLING    PROGRAM. 

C  NfPCIfjNBfXi XB.TEMPV 

C 

C        N  -  THE  SIZE  OF  THE  MATRICES. 

C 

C        MDIM  -  THE  DIMENSIONS  OF  THE  MATRICES  IN  THE  CALLING 

C  PRCGRAM. 

C 

C        NB  -  THE  NUMBER  OF  BANDS  IN  THE  BANDED  MATRIX,  XB. 

C 

C        X  -  TFE  SQUARE  N  BY  N  MATRIX.   DIMENSIONED 

C  (MDIPtMDIM)  IN  THE  CALLING  PROGRAM   (COMPLEX-16) 

C 

C        X3  -  TFE  N  BY  N  MATRIX  WHICH  IS  STORED  IN  BANDED 

C  FORM.   DIMENSIONED  (NB,MDIM)  IN  CALLING  PROGRAM. 

C  (CCMPLEX-16) 

C 

C        TEMPV  -  A  VECTOR  WORKSPACE  WHICH  MUST  BE  PROVIDED  BY 

C  THE  CALLING  PROGRAM.   DIMENSIONED  AT  LEAST  (N). 

C  (CCMPLEX*16) 

C 

C        THE  FOLLOWING  IS  OUTPUT  BY  MULDBM. 

C 

C        X  -  SET  TO  THE  PRODUCT  OF  X  AND  XB  (MDIM, MDIM) 

C  <CCfPLEX*16) 

C 

C      REQUIRED  ROUTINES 

C 

C        NONE 

C 

C 

C 

C 

SUBROUTINE  MULDBM( X, XB , N ,NB, MDIM, TEMPV) 

CC  NP  L  EX*  1  c  X ( MD I M , MD  I M  ) , X B ( N  B , MD I M ) , TE  MP  V ( MD I M ) , T EMP 

NBHM  =  (NB-1J/2 

NBHP  =  (NB  +  D/2 
C 

C  LCCP  CVER  INDEX  I 
C 

CO  100  1  =  1, N 
C 

C  STCRE  COLUMN  I  OF  MATRIX  X  TEMPORARILY 
C 

CC  10  J=1,N 
10  TEMPV( J)  =  X( I  ,J) 
C 

C  FIND  PRODUCTS  FOR  FIRST  NBHM  SPECIAL  CASES,  THAT  IS  WHERE 
C  WHERE  THE  BANDED  MATRIX  DOES  NOT  HAVE  ITS  FULL  WIDTH 


C 


CC    22    J=1,NBHM 
TEMP    =     (CCCODO) 
J J    =    NBHM    +    J 
DO    21    K=1,JJ 
JJJ    =    JJ-K+1 
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c 
c 
c 
c 
c 


c 

c 
c 


21  TEfP    =    TEMP+TEMPV(K)*X8( JJJ,KJ 

22  X(I,J)     =    TEMP 

CGNPUTE    PPCDUCTS    FOR    "REGULAR"    COMBINATIONS    OF    ROWS    AND 
COLUMNS,     THAT     IS,    THOSE    THAT    ARE    NOT    TRUNCATED 
AT    THE    BEGINNING    OR    END    BY    THE    BOUNDARIES 

JF    =    N-NBHM 
CC    32    J=NBFP,JF 
TEMP    =     (CCCODOJ 
CO    31    K=1,N3 
J  J  J    =    NB-K  +  1 

31  TEMP    =    TEMP+TEMPVC  J-NBHP+K)  *X8  (  JJJ  ,J-NBHP+K) 

32  X(IiJ)     =    TEMP 

FIND  PPCDUCTS  FOR  LAST  NBHM  SPECIAL  CASES. 

DC  4  2  J  =  l  ,KBHM 
TEMF  =  (ODCODO) 
JJ  =  NB-J 
CC  41  K=1,JJ 

41  TEMP  =  TEMF+TEMPV(N-JJ+K)*XB(NB-K+1,N-JJ+K) 

42  X(I,N-NBHM+J)  =  TEMP 

IOC  CCNTINUE 

RETURN 
END 
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C SUBROUTINE  MSET 

C 
C 

C  PURPOSE 

c 

C  MSET  SETS  UP  A  MATRIX  FOR  THE  FINITE  DIFFERENCE 

C  PROBLEM  CF  POISEULLE  FLOW  WITH  THE  PROPER  BOUNDARY 

C  CCNDITICNS  FOR  THE  VELOCITY  VECTOR  POTENTIAL  FOR 

C  VISCCLS  FLOW. 

C 

C  LSAGE 

C 

C  CALL  MSET(X,N,MDIM,MODE,CFMAT) 

C 

C  CESCBIPTICN  OF  PARAMETERS 

C 

C  THE  FOLLOWING  MUST  BE  SET  BY  THE  CALLING  PROGRAM 

C  N,MDIN, MODE, CFMAT 

C 

C  N  -  THE  SIZE  OF  THE  MATRIX.   IF  MODE=0,  N  IS  EQUAL 

C  TO  THE  IMUNBEP  OF  POINTS  OF  THE  FINITE  DIFFERENCE 

C  MESH  ACROSS  THE  CHANNEL  NOT  INCLUDING  THE  TWO 

C  BCUNCARIES.   IF  MQDE=1  CR  MCDE=2,  N  IS  EQUAL  TO 

C  CNE-HALF  OF  THE  NUMBER  OF  INNER  POINTS  ACROSS  THE 

C  FULL  CHANNEL.   IN  THIS  CASE,  THE  CHANNEL  MUST  BE 

C  DIVIDED  INTO  AN  EVEN  NUMBER  OF  POINTS  SO  THAT  N 

C  WILL  BE  AN  INTEGER. 

C 

C  MDIM  -  THE  COLUMN  DIMENSION  OF  THE  MATRIX  X.   MDIM 

C  MLST  BE  .GE.  N. 

C 

C  MODE  -  IF  MODE=0,  THE  MATRIX  IS  SET  UP  FOR  THE  FULL 

C  CHANNEL.   IF  MQDE=1  OR  MOUE=2,  THE  MATRIX  IS  SET 

C  UP  FCR  THE  HALF  CHANNEL  AND  THE  BOUNDARY 

C  CONDITIONS  ARE  SET  SUCH  THAT  THE  ODD  OR  EVEN 

C  EIGENFUNCTIONSi  RESPECTIVELY,  ARE  SOLVED  FOR. 

C 

C  CFMAT  -  THE  NAME  OF  A  FUNCTION  SUBPROGRAM  WITH  TWO 

C  PARAMETERS,  K  AND  Y,  INDICATING  WHICH  TERM  OF  THE 

C  FIMTE  DIFFERENCING  IS  DESIRED,  AND  THE  POSITION 

C  RELATIVE  TO  THE  CENTER  OF  THE  CHANNEL.   MUST  BE 

C  DECLARED  EXTERNAL  IN  THE  CALLING  PROGRAM. 

C  (CCMFLEX*16) 

C 

C  THE  FOLLOWING  IS  OUTPUT  BY  MSET 

C 

C  X  -  THE  N  BY  N  MATRIX  INTO  WHICH  THE  COEFFICIENTS 

C  CF  ThE  FINITE  DIFFERENCING  ARE  PUT.   MUST  BE 

C  DIMENSIONED  (MDIM, MDIM)  IN  THE  CALLING  PROGRAM. 

C  (CCMFLEX*16) 

C 

C  NOTES... 

C 

C  THE  EIGENVALUES  AND  VECTORS  OBTAINED  BY  USING  MSET 

C  TWICE,  WITH  MODE=l  AND  MQDE=2,  ARE  THE 

C  SAME  AS  USING  MSET  ONCE  WITH  MODE=0,  BUT  WITH 

C  N  TWICE  AS  LARGE. 

C 

C  OTHER  ROUTINES  NEEDED 

C 

C  NONE  -  EXCEPT  THE  FUNCTION  SUBPROGRAM  NAME  PASSED  IN 

C  THE  PARAMETER  'CFMAT*. 

C 

C 

C 


SUBROUTINE  MSET ( X, N, MD  IM , MODE, CFMAT ) 

RE£L*8  REY,Y,DELY,DFLOAT 

COMPLEX* 16  A,B,G 

CCMMCN  /  CCEFNT  /  A , B , G , REY , DELY 

COMPLEX*  16   CFMAT 

C0MPLEX*16  X(MDIM,MD1M) 
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C  COMPUTE  GRIC  SIZE  FOR  FINITE  DIFFERENCE  MESH  ACROSS 

C  HALF  CHANNEL  CR  FULL  CHANNEL 

C 

DELY    =    2C0/DFL0AT(N+1 ) 

IF(M0DE.£C.l.OR.MODE.EQ.2)     DELY    =    2DO/DFLO AT ( 2*N  + 1 ) 
C 

C    ChECK     IF    MATRIX    DIMENSIONED    LARGE    ENOUGH 
C 

IF(N.LE.MDIM)    GO    TO    1 

WRITE(6,9G00) 
900C    FOPMATPO*    *    *    ERROR   -    ARRAYS    NOT    DIMENSIONED    LARGE1, 
*  *     ENOUGH    *    *    *« ) 

STOP 
C 

C  ZERO  ENTIRE  MATRIX 
C 

1  CO  10  1=1, N 

DC  10  J=1,N 
10  XII , J)  =  (CDO,ODO) 
C 

C  DO  SPECIAL  CASE  AT  DISTANCE  DELY  FROM  CHANNEL  WALL 
C  INCLUDING  BCUNDARY  CONDITIONS 
C 

Y  =  1D0-DELY 

X(l,l)  =  CFMAT(3,Y)+CFMAT( 1,Y) 

X(l,2)  =  CFMAT(4,Y) 

X(i,3)  =  CFMAT(5,Y) 
C 

C    DO    SPECIAL    CASE    AT    DISTANCE    2*DELY    FROM    CHANNEL    WALL 
C    INCLUDING    ECUNCARY    CONDITIONS 
C 

Y  =     1DC-2C0*DELY 
X(2,l)     =    CFMAT(2,Y) 
X(2,2j     =    CFMAT(3,Y) 
X(2,3)    =    CFMAT(4,Y) 
X(2,4)     =    CFMAT15.Y) 

C 

C  DO  ALL  REGULAR  POINTS  IN  BETWEEN,  THAT  IS,  THOSE  VALUES 

C  OF  Y  FCR  WHICH  ALL  5  FINITE  DIFFERENCE  GRID  POINTS  ARE 

C  INTERICR  TC  THE  CHANNEL 

C 

IL  =  N-2 

CO  20  1  =  3, IL 

K  =  1-3 

Y  =  1CC-CELY*DFLCAT( I) 
DO  20  J=l,5 

20  X( I,K+J)  =  CFMAT( J,YJ 
C 

C  FINALLY  DG  THE  TWO  SPECIAL  CASES  WHICH  OCCUR  EITHER  AT 
C  THE  CENTER  OF  THE  CHANNEL  OR  AT  THE  OTHER  WALL,  DEPENDING 
C  CN  THE  VALUE  OF  MODE.   BOUNDARY  CONDITIONS  ARE  SET  UP 
C  DEFENDING  CN  MODE 
C 

Y  =  1CC-CELY*DFL0AT(N-1) 
X(N-l,N-3)  =  CFMATtl ,Y) 
X(N-l,N-2)  =  CFMAT(2,Y) 
X(N-1,N-1)  =  CFMAT(3,Y) 
XIN-1,N)  =  CFMAT(4,Y) 

IF(MCDE.EQ.l)     X(N-1,N)     =    CFMAT ( 4 ,Y ) -CFMA T( 5, Y) 
IFIMCDE.EQ.2)     X(N-1,N)     =    CFMAT (4 ,Y ) +CFMATI 5, Y ) 
C 

Y  =     1CC-DELY*DFL0AT( N) 
X(N,N-2)     =    CFMAT( 1 ,Y) 
X(NiN-l)     =    CFMAT(2,Y) 

IF(MCDE.EQ.l)  X(N,N-1)  =  CFMAT ( 2 ,Y } -CFMA T( 5, Y) 
IF(MQD£.EQ.2)  X(N,N-1)  =  CFMAT ( 2 ,Y ) +CFMAT( 5, Y ) 
X(N,N)     =    CFMAT(3,YJ+CFMAT(5,Y) 

IF(MODE.EQ.l)     X(N,N)     =    C FM AT ( 3 , Y ) -CFMAT( 4, Y ) 
IF(MCDE.EG.2)    X(N,N)     =    CFMAT ( 3 , Y )+CFMAT { 4, Y J 
C 

RETURN 
END 
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C SUBROUTINE    BMSET 

C 
C 

C      PURPOSE 

c 

C        THE  PURPOSE  OF  BMSET  IS  EXACTLY  THAT  OF  MSET  EXCEPT 

C        THAT  BMSET  TAKES  ADVANTAGE  OF  THE  BANDED  NATURE  OF 

C        THE  FINITE  DIFFERENCE  MATRICES  TO  CONSERVE  SPACE. 

C 

C      LSAGE 

C 

C        CALL  BMSET(X,N,MDIM, MODE, CFMAT) 

C 

C      DESCRIPTION  OF  PARAMETERS 

C 

C        THE  PARAMETERS  FOR  BMSET  ARE  THE  SAME  AS  THOSE  FOR 

C        MSET  kITH  THE  EXCEPTION  THAT  THE  MATRIX  X  MUST  BE 

C        DIMENSIONED  (5,MD1M)  IN  THE  CALLING  PROGRAM. 

C 

C      NOTE:   THE  PROCEDURE  IS  IDENTICAL  TO  THAT  OF  MSET. 

C        COMMENTS  HAVE  THEREFORE  NOT  BEEN  INCLUDED  IN  BMSET. 

C 

C 

c 

c 

SUERCUTINE    BMS ET (X ,N , MD I M, MODE ,C FMAT ) 

PEAL*8    REY,Y,DELY,DFLOAT 

CCMPLEX-16       CFMAT 

CCMPLEX*16    A,B,G 

CCMPLEX*16    X(5,MDIM) 

COMMCN    /    CCEFNT    /    A, B , G ,REY, DELY 

CELY    =    2DO/DFLCAT(N+l) 

2DO/DFLOAT(2*N+l) 

9000    FORMAT  (  «' 0-    *    *    ERROR    -    ARRAYS    NOT    DIMENSIONED    LARGE*, 


1C 


IF(MCDE. 

EC. 1. OR. J 

^10DE 

.EQ.2) 

DELY    = 

I  F  (  N  .  L  E  . 

f'.DIM)    GO 

TO 

1 

kRITE(6, 

9000) 

FORM  AT ( « 

0^    *    *    El 

*ROR 

-    ARRAY 

■     ENCUGH    *    * 

*•  ) 

STOP 

CO    10    1= 

If  5 

CO    10    J= 

1,.VDIM 

X(I , J)     = 

(CUO,ODO) 

Y    =     100- 

DELY 

X(3,l)    = 

CFMAT (3 

,Y)  + 

CFMAT( 1 

,Y) 

x(4,n  = 

CFMAT (4 

,Y) 

X(5,l)    = 

CFMAT (5 

,Y) 

Y    =     1CC- 

2D0*DELY 

xta,2)   = 

CFMAT  (2 

,Y) 

X(3,2)     = 

CFMAT (3 

,Y) 

X14,2)     = 

CFMAT (4 

,Y) 

X(5,2)     = 

CFMAT (5 

,Y) 

IL    =    N-2 

• 

CO  20  1  =  3, IL 

Y  =  1CC-DELY*DFL0AT( I) 
CO  20  J=l,5 

2C  X( J, I)    =  CFMATl J,Y) 

Y  =    1DO-DELY#DFLOAT(N-1) 
X(1,N-1)  =    CFMAT(1 ,Y) 
X(2,N-1J          =    CFMAT(2,Y) 
X(3,N-1)          =    CFMAT(3,Y) 
X(4tN-l)     =    CFMAT(4,Y) 

IF(MCDE.EQ.l)     X(4,N-1)     =    CFMAT ( 4 , Y )-CFMAT( 5, Y) 
IFCM0DE.EG.2J     X(4,N-1)    =    C FMAT ( 4 ,Y ) +CFMAT< 5, Y) 

Y  =     1CC-DELY*DFL0AT(NJ 
X(  1,N)  =    CFMATt 1,Y) 
X(2,N)  =    CFMAT(2,Y) 

IF(MCDE.EC.l)     X(2,N)     =    CFM AT ( 2 , Y )-CFMAT ( 5, Y ) 

IF(M0D£.EQ.2)  X(2,N)  =  CFMAT( 2 , Y ) +CFMAT( 5, Y) 

X(3,N)  =  CFMAT(3,Y )+CFMAT( 5,Y) 

IF(MODE.EQ.l)  X(3,N)  =  CFM AT < 3 , Y ) -CFMAT ( 4, Y ) 

IF(MCDE.EU.2J  X(3,N)  =  C FMAT ( 3 , Y ) +  CF MAI i 4, Y ) 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c. 

c 


SUBROUTINE  DSPLIT 


PURPOSE 


DSPLIT  TAKES  A  MATRIX  OF  CCMPLEX-16  NU 
SPLITS  IT  INTO  TWO  MATRICES,  ONE  CONTA 
PART  CF  THE  ORIGINAL  MATRIX,  AND  ONE  C 
IMAGINARY  PART. 

USAGE 

CALL  CSPLIT(N,MDIM,A,AREAL,AIMAG) 

DESCRIPTION  OF  PARAMETERS 

N  -  THE  SIZE  OF  THE  MATRIX  A,  AN  N  BY 
MATRIX. 


MBERS  AND 
INING  THE  REAL 
ONTAINING  THE 


N  SQUARE 


MDIM  -  THE  COLUMN  DIMENSION  OF  MATRIX  A 


A  -  THE  INPUT  MATRIX.  MUST  BE  DIMENSI 
AT  LEAST  N  IN  THE  CALLING  PROGPAM   ( 

AREAL, AIPAG  -  THE  OUTPUT  MATRICES  CONT 
REAL  AND  IMAGINARY  PARTS,  RESPECTIVE 
MATRIX  A.  MUST  BE  DIMENSIONED  (MDIM 
CALLING  PROGRAM. 

NOTES. .. 

MATRIX    A    AND    MATRIX    AREAL    MAY    OVERLAP 
DIMENSIONED    IN    THE    CALLING    PROGRAM    A 


ON  ED    MDIM    BY 
C0MPLEX*16) 

AINING    THE 

LY,     OF 

t MDIM J     IN    THE 


IF    THEY    ARE 
S    FOLLOWS... 


CCMPLEX*16    A (MDIM, MDIM) 

REAL*8    AREAL (MDIM, MDIM) , A  I  MAG ( MDI M , MO I M ) 

EQUIVALENCE (A(l ,1) ,AREAL(  1,1)) 


CTFER  ROUTINES  NEEDED 
NONE 


SUBROUTINE  DSP  LI T( N, MD IM , A , AR, A I ) 

REAL*8  A(2t MDIM, MDIM) , ARC MDIM, MDIM J , AI (MDIM, MDIM) 

DC  1  J=1,N 
CO  1  1=1,  N 
AR(I  ,JJ  =  Ad,  I, J) 
1  AKI  ,J)  =  A(2,I,  J) 

RETURN 
END 
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C SUBROUTINE  DEIGEO 

C 
C 

C      FUFFOSE 
C 

C  DEIGEC    SOLVES    THE    LINEARIZED    N AV IER-ST OKES    EQUATION 

C  FOR    POISEULLE    FLOW.       THE     INPUTS    TO    DEIGEO    ARE    THE 

C  WAVE    NCS.i     ALPHA    AND    BETA,     AND    THE    REYNOLDS    NO. 

C  DEIGEC    OUTPUTS    THE    EIGENVALUES    FOR   GAMMA. 

C 

C  USAGE 

C 

C  CALL    CEIGEOt ALPHA, BETA, REYNC,N ,MDI M, WREVEN , WIEVEN, 

C  WRCDD,WIODD,CDM,DM,BM,WV) 

C 

C      CESCFIPTICN  OF  PARAMETERS 

C 

C        THE  FOLLOWING  MUST  BE  SET  BY  THE  CALLING  PROGRAM... 

C  ALPHA, BETA, REYNO.N, MDIM 

C 

C        ALPHA  -  THE  PERTURBATION  WAVE  NUMBER  IN  THE  FLOW 

C  DIRECTION  (X).   (CCMPLEX*16) 

C 

C        BETA  -  THE  PERTURBATION  WAVE  NUMBER  IN  THE  DIRECTION 

C  PERPENDICULAR  TO  THE  FLOW  BUT  PARALLEL  TO  THE 

C  PLATES  (ZJ.   (C0MPLEX*16) 

C 

C  REYNC    -    THE    REYNOLDS    NUMBER    (REAL*8) 

C 

C        N  -  THE  SIZE  OF  THE  MATRICES  WHICH  IS  EQUAL  TO 

C  (NC-D/2  WHERE  ND  IS  THE  NUMBER  OF  DIVISIONS 

C  ACRCSS  THE  CHANNEL.   (NOTE...  DEIGEO  SOLVES  THE 

C  PROBLEM  ACFOSS  THE  HALF  CHANNEL  TWICE  -  ONCE  FOR 

C  THE  EIGENVALUES  CORRESPONDING  TO  THE  EVEN 

C  EIGENFUNCTIONS  AND  ONCE  FCR  THOSE  CORRESPONDING 

C  TO  THE  ODD  EIGENFUNCTIONS. 

C 

C        MDIM  -  ThE  COLUMN  DIMENSION  OF  THE  MATRICES  WHICH 

C  DEIGEO  USES.   MDIM  MUST  BE  . GE .  N 

C 

C        THE  FOLLOWING  ARE  OUTPUT  BY  DEIGEO 

C  WREVEN, WIEVEN, WRCDD.WIOOD 

C 

C        WREVEN, WIEVEN  -  THE  REAL  AND  IMAGINARY  PARTS  OF  THE 

C  EIGENVALUES  CORRESPONDING  TC  THE  EVEN 

C  EIGENFUNCTIONS.   DIMENSIONED  TO  AT  LEAST  N. 

C  (REAL*8) 

C 

C  WRCDD,WIODD    -    THE    REAL    AND    IMAGINARY    PARTS    OF    THE 

C  EIGENVALUES    CORRESPONDING    TO    THE    ODD 

C  EIGENFUNCTIONS.       DIMENSIONED    TO    AT    LEAST    N. 

C  (REAL*8) 

C 

C  THE    FCLLCWING    MATRICES    MUST    BE    INPUT    TO    DEIGEO    AS 

C  hCPKSPACE. 

C 

C  CDM(MDIM,MDIM)       (CCMPLEX*16) 

C  CM(MCIM,MDIM)       (REAL -8) 

C  BM(5,MDIM)        (C0MPLEX*16) 

C  WV(MDIM)        (CCMPLEX-16) 

C 

C  NOTES... 

C 

C  THE    MATRICES    CAN    BE    OVERLAPPED    TO    CONSERVE    SPACE, 

C  FCR    EXAMPLE,    FOR    N    =    60... 

C 

C  CCMPLEX*16    CDM(60,30,3) 

C  CCMPLEX*16    DM(60,30), WV(60) ,BM(5,60) 

C  EQUIVALENCE     ( DM ( 1 , 1 ) ,CDM(  1 ,  1 , 3 )  ) , 

C  (13,111,1  ),CDM  1,  1,31  ), 

C  WVEC(l)  ,CDM(1  ,6,3)  ) 

C 
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C        NOTE  THAT  IT  IS  ONLY  THE  ACTUAL  SIZE  OF  THESE 

C        WORKSPACES  THAT  IS  IMPORTANT,  NOT  THEIR  TYPE. 

C 

C      CTHER  FCLTINES  NEEDED 

C 

C        THE  FCLLCWING  ARE  CALLEC  BY  DEIGEO 

C 

C  CHM1E1,CHM2E2,MSET,CDMTIN,BMSET,MULDBM,DSPLIT, 

C  EFESSCELRH1C 

C 

C 

C 

SUBROUTINE  DEI GE C( AL PH A, BETA , REYNO ,N ,MDI M, 
*     WR EVEN, WIEVEN.WR ODD, WIODD, CDM, DM,BM,WV) 

IMPLICIT  CCMPLEX*16( A-H, O-Z) 

DIMENSION  IVEC(IOO) 

REAL*8  ksREVEN(l)  tWIEVEN(  1)  ,WRODD(l  )  ,  WIODD(  1) 

PEAL*8  CCMC1 ),DM<1 ),BM(1 ),WV(1) 

FEAL*8  PEY,DELY, REYNO 

COMMON  /  CCEFNT  /  A  ,  B , G ,RE Y,D£LY 

EXTERNAL  ChMlEl , CHM2E1 
C 

C  THIS  SUBROUTINE  SOLVES  THE  EQUATION   YV  =  GXV   WHERE 
C  X  ANC  Y  ARE  MATRICES,  V  IS  THE  EIGENVECTOR  AND  G  IS  THE 
C  EIGENVALUE.   TEE  EIGENVALUES  ARE  DETERMINED  AND  PASSED 
C  BACK  TO  TEE  CALLING  PROGRAM  IN  WRODD,  WIODD,  WREVEN  AND 
C  WIEVEN. 
C 

A  =  ALPHA 

E  =  BETA 

REY  =  PEYNC 
C 

C  SET  LP  MATRIX  X  FOR  ODD  EIGENVECTORS. 
C 

CALL  MSET(CDM,N,MDIM,1 ,CHM2E1) 
C 

C  INVEPT  MATRIX  X. 
C 

CALL  CCMTIN(N,CDM,MDIM,DETERM) 
C 

C  SET  UP  MATRIX  Y  IN  BAND  STORAGE  MODE  FOR  ODD  EIGENVECTORS. 
C 

CALL  BNSET(BM,N,MDIM,1,CHM1E1) 
C 

C  MULTIPLY  MATRIX  Y  BY  THE  INVERSE  OF  MATRIX  X  TO  CONVERT 
C  TO  THE  STANCARC  EIGENVALUE  PROBLEM  WHICH  HAS  THE  FORM 
C  (Z-GJV  =  0   WHERE   Z  =  ( Y) ( I NVE RSE ( X ) ) . 
C 

CALL  MUL C BM ( CDM , BM , N , 5 , MD I M , WV  ) 
C 

C  SPLIT  MATRIX  INTO  REAL  AND  IMAGINARY  PARTS  AND  CALL 
C  THE  SLSRCUTINES  TO  FIND  THE  EIGENVALUES. 
C 

CALL  DSPLIT(N,MDIM ,CDM,CDM,DM) 

CALL  EHESSC( CDM, DM,1 ,N ,N ,MDI M, I VEC  ) 

CALL  ELPH1C(  CDM, DM,1 ,N,M,MDIM, WPODD, WIODD, IN  ERR, IER) 

IFUNERR.NE.O)  WRI  TE  (  6  ,9000  J  INERR,IER 
90CC  FCRMATCOERROR  NUMBER', 17,'   ON  EI  GENV  ALUE  •  ,  I  7,  ///  ) 
C 

C  REPEAT  TEE  SOLUTION  FOR  EIGENVALUES  FOR  THE  EVEN 
C  EIGENVECTORS 
C 

CALL  MSET(CDM,N,MDIM,2,CHM2E1) 

CALL  C CMT I N( N , CD  f , MD I M , D ET  E RM ) 

CALL  BMSET(BM,N,MDIM,2,CHM1E1) 

CALL  MULLBMICDM,BM,N,5,MDIN,WV) 

CALL  D  S  P  L I T ( N , MD I M , C  DM , C  DM , D  M ) 

CALL  EhESSC( CDM, DM,1 ,M,N,MDIM, IV EC) 

CALL  ELRH1C(CDM, DM, 1 ,N,N,M DIM, WREVEN t WIEVEN, INERR, IER) 

IF( INEHK.NE.O)  WRITE(6,9000)  INERR, IER 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c, 

c 


c 
c 
c 


c 
c 
c 


SUBROUTINE    CDMTIN    (CATEGORY    F-l ) 
PURPOSE 

INVERT    A    C0MPLEX*16    MATRIX 
USAGE 

CALL    CDMTIN(N, A, NDIM, DETERM) 

DESCRIPTION    OF    PARAMETERS 

N  -    ORDER    OF    COMPLEX*l6    MATRIX    TO    BE    INVERTED 

(INTEGER)     MAXIMUM    «N«     IS    100 

A  -    C0MPLEX*16    INPUT    MATRIX    (DESTROYED).       THE 

INVERSE    OF     «A«     IS    RETURNED    IN    ITS    PLACE 

NDIM         -    THE    SIZE    TO    WHICH    'A'     IS    DIMENSIONED 

(ROW    DIMENSION    OF    'A*     ACTUALLY    APPEARING 
IN    THE    DIMENSION    SFATMENT    OF    USER'S 
CALL  ING    PROGRAM) 

DETERM    -    C0MPLEX*16    VALUE    OF    DETERMINANT    OF     «A» 
RETURNED    BY    CDMTIN. 

REMARKS 

MATRIX    «A«     MUST    BE    A    C0MPLEX*8    GENERAL    MATRIX 

IF    MATRIX    «A«     IS    SINGULAR    THAT    MESSAGE    IS    PRINTED 

»N«     MUST    BE    .LE.    NDIM 

SUBROUTINES    AND    FUNCTIONS    REQUIRED 

ONLY    BUILT-IN    FORTRAN    FUNCTIONS 

METHOD 

GAUSSIAN  ELIMINATION  WITH  COLUMN  PIVOTING  IS  USED. 
THE  DETERMINANT  IS  ALSO  CALCULATED.  A  DETERMINANT 

OF  ZERO  INDICATES  THAT  MATRIX  'A'  IS 

SINGULAR. 


SUEPCUTINE  CDMTIN  ( N , A ,MDI M, DETE RM ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMPLEX* 16  A (NDIM, NDIM), PIVOT (100) ,AMAX, T, SWAP, 

4t  P  C  T  C  D  JJ   II 

INTEGERS  IPIVOT(IOO),  INDEX(  100,2) 
REAL*8  TEMP,ALPHA(10C) 

INITIALIZATION 

DETEFM  =  (IDOtODO) 

DC  20  J=1,N 

ALFHA(  J)  =  ODO 

DC  10  1=1, N 
10  ALPHA(J)=ALPHA(J)+A(J, I)*DCCNJG(A( J, I) ) 

ALPHA ( J) =DSQRT( ALPHA ( J) ) 
20  IFIVOT(J)=C 

DC  600  1  =  1, N 

SEARCH    FOR    PIVOT    ELEMENT 


ANAX    =     (ODO, ODO) 
CO    1C5    J=l,.M 

IF     (  IPIVCT(J)-l)    60, 105,60 
60    CO    100    K=1,N 

IF    ( IPIVOT(K)-l)    80,100,740 
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8C    TEMP=ANAX*DCONJG(AMAXJ-A( J,K)*DCONJG(A(J,K)) 

IF(TEMP)     £5,     85,     100 
£5     IRC^=J 
ICCLUN=K 
AMAX=A ( J,K) 
ICC    CCMINLE 
105    CONTINUE 

IPIVCT  (ICGLUM)=IPIVOT( IC0LUM)+1 
C 

C  INTERCHANGE    ROWS    TO    PUT    PIVOT    ELEMENT    ON    DIAGONAL 

C 

IF( IRCW-ICCLUM)     140,     260,     140 
140    CETERM=DETEP.M 
DC    200    L  =  1,N 
SWAF=A(IFCW,L) 
A( IKOW,L)=A( ICOLUM,L) 
20C    A(  ICCLLN,L)=SWAP 
SWAP=ALPHA(IROW) 
ALPHA!  IRCW)=ALPHA( ICOLUM) 
ALPHAC  ICCLUM)=SWAP 
260    INCEXi  I,1)  =  IR0W 

IND5X(  1, 2)=IC0LUM 
PIVGTC I )=A( ICOLUM, ICOLUMJ 
L=PIVOT<  I) 

TcMP=PIVCT(I)*DCCNJG(PIVOT( I )) 
IF(TEMP)    330,     720,    330 
C 

C       DIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 
C 

330  A( ICOLUM, ICOLUM)  =  (1D0,0D0) 
CC  350  L=1,N 
L=FIVOT(  I) 
350  A( ICOLUM,L)=A( ICOLUM,L)/U 
C 

C       REDUCE  NCN-PIVOT  ROUS 
C 

38C  CC  550  L1=1,N 

IF(Ll-ICCLUM)  400,  550,  400 
400  T=A(L1,ICCLUM) 

A(L1,ICCLUMJ  =  (0D0,0D0) 
CC  450  L=1,N 
L  =  A(  ICCLU^,L) 
450  A(L1,L)=A(L1,L)-U*T 
55C  CCMINLE 
6CC  CONTINUE 
C 
C       INTERCHANGE  COLUMNS 


C 


62C  DC  71C  1  =  1, N 

L=N+1-I 

IF( INDEX(L,1 )-INDEX(L, 2) )  630,  710,  630 
63C  JR0W=INDEX(L,1) 

JCCLUM=INDEX(L|2J 

DO  7C5  K=1,N 

SWAP=A(K,JPOW) 

A(K, JPCW)=A( K, JCCLUM) 

A(K, JCCLUMJ^SWAP 
7C5  CONTINUE 
710  CCNTINUt 

RETURN 
720  WRITE(6,73C) 

730  FCRMAT(20H   MATRIX  IS  SINGULAR) 
74C  FETUFN 

END 
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C SUBROUTINE  DFIT2 

C 
C 
C      PURPOSE 

c 

C  THIS  SUBROUTINE  TAKES  3  POINTS  (XI, Y),  (X2tY2), 

C  AND  (X2,Y3),  AND  USES  A  PARABCLIC  FIT  TO  FIND  THE 

C  EXTREMAL  VALUE  OF  Y  AND  THE  CORRESPONDING  VALUE  FOP 

C  X  WHICH  ARE  RETURNED  IN  XM  AND  YM.   IT  ALSO  TAKES  AN 

C  INPUT  VALUE,  YVAL,  AND  COMPUTES  THE  TWO  SOLUTIONS 

C  FOR  X  TO  THE  PARABOLIC  INTERPOLATION,  TAKES  THE 

C  SOLUTICN  NEAREST  TO  XI,  AND  RETURNS  IT  IN  XNEXT. 

C  IF  TFEPE  IS  NO  REAL  SOLUTIGN  TO  THE  EQUATION,  THEN 

C  XNEXT  IS  SET  TO  O.O. 

C 

C  USAGE 

C 

C  CALL  CFIT2(X1,X2,X3,Y1,Y2,Y3, YVAL, XNEXT, XM,YM) 

C 

C  DESCRIPTION  OF  PARAMETERS 

C 

C  THE  FOLLOWING  MUST  BE  SET  BY  THE  CALLING  PROGRAM 

C  X1,X2,X3,Y1,Y2,Y3,YVAL 

C 

C  X1,X2,X3,  -  THE  X  VALUES  OF  3  POINTS  (REAL*8) 

C  Yi,Y2,Y3  -  THE  Y  VALUES  OF  3  POINTS  (REAL*8) 

C  YVAL  -  A  VALUE  OF  Y  FOR  WHICH  THE  CORRESPONDING 

C  X  VALUE  IS  DESIRED  (REALMS) 

C 

C  THE  FOLLOWING  ARE  SET  BY  DFIT2 

C  XNEXT, XMiYM 

C 

C  XNEXT  -  SET  TO  THE  VALUE  OF  X  NEAREST  XI 

C  CORRESPONDING  TO  YVAL  I REAL*£) 

C  XMtYM  -  SET  TO  THE  MIN  Y  VALUE  FROM  THE  SECOND 

C  GRCER  INTERPOLATION,  WITH  THE  CORRESPONDING  X 

C  VALUE  (REAL*8) 

C 

c 

c 

c 

SUBROUTINE  DFIT2 (X 1 , X2 ,X3, Y 1 , Y2 , Y3 , YVAL, XNEXT,XM,YM) 

IMPLICIT  REAL*8(A-Z) 
C 

A  =  ( (Y1-Y2)/(X1-X2)-(Y1-Y3)/(X1-X3) ) 
*     /( (Xl**2-X2-*2)/(Xl-X2)-(Xl**2-X3**2)/(Xl-X3) ) 

B  =  (Y1-Y2-A*(X1**2-X2**2))/(X1-X2) 

C  =  Yl-A*Xl**2-d*Xl 

XM  =  -B/(2D0*A) 

YM  =  A*XM**2+S*XM*C 

D  =  B**2-4D0*A*(C-YVAL) 

IF(D.GE.CDO)  GO  TO  10 

XNEXT  =  CDO 

RETURN 
C 

1C  XN1  =  (-B+DSQRT(D) )/(2D0*A) 

XN2  =  (-B-DSQRT(D)  }/(2D0*A) 

XNEXT  =  XN1 

IF(DABS(XN2-XI) .LT.DABS(XNl-Xl) )  XNEXT  =  XN2 

RETURN 

END 
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c . 

C  PFCGFAM  #R1A 

C 

C  PROGRAM  TO  PRINT  EIGENVALUES 

C  FOR  THE  3-D  CYLINDRICAL  FLOW  PROBLEM 

C  N  =  0       FOR  G 

C 

C 

c 

C        INPUT  VALUES  ARE  N,  REYr  AR,  AND  AI. 

C  (USING  NAMELIST  'LIST* ) 

C 

c 

IMPLICIT  REAL*8< A-H, 0- Z ) 

CCMPLEX*16  A,B,DETERM,G 

CC MP  LEX- 16  XMATC31  ,31)  ,YMAT(31,31)  ,WV(31) 

CCMPLEX*16  CGM1E2,CGM2E2 

FEAL*8  GR(31) ,GI(31) 

REAL*4  GR4(31) ,GI4(31) 

CIMENSICN  IVEC(IOO) 

CCMMCN  /  CCEFNT  /  A, 8, G,REY, DELR 

NAMELIST  /  LIST  /  N,REY,AR,AI 

EXTERNAL  CGM1E2,CGM2E2 
C 

NCIM  =  31 

N  =  30 

REY  =  6C00DO 

AR  =  OCO 

AI  =  1D0 

B  =  (ODOfODO) 
C 

1  READ(5,LIST,END=100) 

A  =  DCMPLX(AR,AI ) 

WR1T5(6,9C00) 
900C  FCRMATf  1«  ) 

hRIT£< 6,9001)  N,REY,A 

9001  FCRMAT(e  CCMPONENT  G  IN  EQUATION  ZS/.'ON  =',14,/, 

*  •  FEY^r-lO^SX,' ALPHA  =  •  ,  2F1  2.  7  ,  8X  ,  •  NBI  =  0s) 
C 

NP1  =  N+l 

CALL  SET2(XMAT,NP1 ,1  ,  MDI M, CGM2E2 , 1 ,1) 
CALL  CCMTIN(N,XMAT,MDIM,DETERM) 
CALL  SET2(YMAT,NP1,1 , MDI M, CGM1 E2 , 1 ,1) 
CALL  MU  L  M ( XM A  T , Y  MA  T , N , MD I M , U V ) 
CALL  DSPLIT(N,MDIM,XMAT,XMAT,YMAT) 
CALL  EFESSC(XMAT,YMAT,  1 ,  N  ,  N  ,  MDI  M  ,  I  VEC) 
CALL  ELRHiC(XMAT,YMAT,l,N,N,MDIM,GP,GI,INERR,IER) 
IF(iNERR.NE.O)  WRI TE (6  ,  9010 )  INERR,IER 
9010  FORMAT(«0*  *  *  ERROR  NUMBER', 17,'  ON  EIGENVALUE', 

*  17 »•  *  *  *«,///) 
C 

WRITE(6,9CC2) 

9002  FCRMAT(///,16X,'GAMHA    REAL ', 5X ,' GAMMA    IMAGM 
CC    10     1=1, N 

GR4(I)     =    SNGL(GRd)) 
GI4(  I)    =    SNGHGI  (I  )) 
10    kRITE(6,9003)     I ,  GR  (  I  ) ,  GI  ( I  ) 

9003  FCFMAT(«0« fI10flF2D15.4) 


C 


WRITE(6,9C00) 

CALL  FLCTP(GR4,GI4,N,0) 

V»RITE(6,9C01  )  N,REY,A 

GO  TO  1 
IOC  WRITE(6,90CO) 
STOP 
END 
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c , 

C  PFCGFAM  #R16 

C 

C  PROGRAM  TO  PRINT  EIGENVALUES 

C  FOR  THE  3-D  CYLINDRICAL  FLOW  PROBLEM 

C  N  =  0       FOR  H 

C 

c , 

c 

C        INPUT  VALUES  ARE  N,  REY,  AR,  AND  AI. 

C  USING  NAMELIST  'LIST' ) 

C 

C o , 

IMPLICIT  REAL*8( A-H,C-Z) 

CCMPLEX*I6  A,S,DETERM,G 

COMPLEX*  16  XMAT(31  ,31  )  ,YMAT(31,31)  ,WV(31) 

CCMPLEX*16  CHM1E1.CHM2E1 

REAL-8  GR(31),GI(31) 

REAL*4  GR4<31) ,G14(31) 

DIMENSION  IVEC(IOO) 

COMMON  /  CCEFNT  /  A  ,  B , G,REY, DELR 

NAMELIST  /  LIST  /  N,REY,AR,AI 

EXTEFNAL  CHM1 El , CHM2E1 
C 

MDIM  =  31 

N  =  30 

REY  =  6C00D0 

AR  =  OCO 

AI  =  ICO 

B  =  (OCO, OCO) 
C 

1  PEAD(5,LIST,END=100) 

A  =  DOFLX(AR,AI) 

WRITE(6,9C00) 

900 0  FORMAT («  1* } 
WRIT£(6,S001 )  N,REY,A 

9001  FCRM^TC  COMPONENT  H  IN  EQUATION  1 

*  •  REY=',F10.2,8X,' ALPHA  =• , 2F12. 7, 8X  ,  • N3I  =  0* 
C 

NP1  =  N+l 

CALL  SET2<XHAT,NPI ,1 , MDI M, CHM2E 1  ,1  ,1) 
CALL  CCMT INC  N,XMAT ,MCIN, CETERM) 
CALL  SET2(YMAT,NP1 ,1 , MDI M, CHM1E1 , 1  ,1 ) 
YMAT(N,N)  =  YMAT(N,N)  +  CHM1 El ( 5 ,DELR) 
CALL  MULMC  XMAT ,YMAT , N ,MDIM , WV ) 
CALL  DSPLIT(N,MDIM,XMAT,XMAT,YMAT) 
CALL  EFESSC(XMAT,YMAT, 1 , N, N ,MDI M , I VEC) 
CALL  ELRH1C(XMAT,YMAT, 1 , N, N , MD I M ,GR, GI ,INERR,IER) 
IF< INERP.NE.O)  WBITE(6,9010)  INERR,IER 
9010  FCRMATPO*  *  #  ERROR  NUMBER1, 17,  ■  ON  EIGENVALUE1, 

*  17,  «  *  *  *• ,///) 
C 

WRITE (6, 9002) 

9002  FCRMAT(///,16X,« GAMMA  REAL •  , 5X, • GAMMA  IMAG* ) 
DC  10  1=1, N 

GR4< I)  =  SNGL(GR(I )) 
GI4( I)  =  SNGLiGI (I )) 
10  WRITE{6,9003)  I , GR ( I ) , GI ( I ) 
90C3  FCRMATCO*  ,  1 10  ,  1P2D15  .4) 

WFITE(6,9000) 

CALL  FLOTP(GR4,GI4,N,0) 

WRITE(6,9001 )  N,REY,A 

GC  TO  1 
100  ^RITE(6,9000) 

STOP 
END 


C 
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c 

C  PROGRAM  ^RZ 

C  PROGRAM  TO  PRINT  EIGENVALUES 

C  FOR  THE  3-D  CYLINDRICAL  FLOW  PROBLEM 

C  N  =  1 

C 

C 

c 

IMPLICIT  REAL*8( A-H,C-ZJ 

CCMPLEX*16  AiB,DET£RM,G 

COMPLEX* 16  XMAT( 60,60) ,YMAT(60,60) ,WV(60) 

CQMPLEX*16  CFMlEl,CFM2ElfCFMl£2»CFM2E2,CGMlEl,CGM2El, 

*  CGMlE2,CGM2E2,CHMlEl,CHM2El,CriMlE2,CriM2E2 
REALMS    GRUO)  ,GI(6G) 

PEAi_*4    GR4(60)  ,GI4(60) 

CIMENSICN  IVECdOOJ 

CCMMCN  /  CCEFNT  /  A , 8 , G, REY ,DELR 

NAMELIST  /  LIST  /  N,REY,AR,AI 

EXTERNAL    CFM1E1 , CFM2 El ,CFM1 E2, CFM2E2 ,C GM IE  1 ,CGM2E1 , 

*  CGM1E2,CGM2E2,CHM1E1 ,CHM2E1 , CHM1E2,CHM2E2 
C 

MDIM  =  6C 
N  =  30 

FEY    =    6C00D0 
AR    =    OCO 
AI    =     ICO 
B    =    (0CC1D0J 
C 

1    REAC(5,LIST,END=100) 
A    =    CCMPLX(AR,AI ) 
WRITE(6,9000) 
90CC    FCFMAT (• 1' ) 

WRITE(6,9001)    N,REY,A 

9001  FCRMATC  N  =',14,/,'  REY  =  •  ,  F10  .  2  ,  8X  ,  •  ALPHA  =',2F12.7, 

*  8X,'NBI  =  1',/,'  FOR  VELGCITY  VECTOR  POTENTIAL  «, 

*  «H  SET  EQUAL  TC  ZERO' ) 
C 

MSIZE  =  2*{N-1) 

CALL    SET2<XMAT,N,1 ,MDI M, CFM2E1 , 1,1) 

CALL    SET2<XMAT,N,1,MDIN,CGM2E1,1,N) 

CALL    SET2(XMAT,N,1 ,MDI M , CFM2 E2, N , 1 ) 

CALL    SET2(XMAT,N,1 ,MCI^,CGM2E2,N,N) 
C 

XMAT(N-1,N-1 )     =    XMAT(N-1,N-1 )     +    CFM2 El ( 5 , DELR ) 

Xf.AT(MSIZE,N-l)     =    XMAT(MSIZE,N-1)     *    CFM2  E2  C  5  ,  DELR  ) 
C 

CALL    CDMTIN(MSIZE,XMAT,MDIM,DETERM) 
C 

CALL     SET2(YMAT,N,1 ,MOI P, CFM1E1 , 1,1  J 

CALL     SET2(YMAT,N,1 ,MDI M, CGM 1E1 , 1,N) 

CALL    SET2(YMAT,N,1 ,MDI M , CFM1 E2 , N, 1 ) 

CALL    SET2(YMAT,N,1,MDI^,CGM1E2,N,N) 
C 

YMAT(N-1,N-1 )     =    YMAT(N-1 ,N-1 )    +    CFM1 El ( 5 ,DELR) 

YMAT(MSIZE,N-1)     =    YMAT( MSI ZE , N-l )     +    CFM1 E2 (5 , DELRJ 
C 

CALL    MULMXMAT,YMAT,  MSI  ZE,  MDIM,  WV) 

CALL    DSPLIT(MSIZE,MDIM,XMAT,XMAT,YMAT) 

CALL    EhESSC(XMAT,YMAT, 1 , MSI ZE ,MSIZE , MDI M , I VEC) 

CALL    ELRH1C(XMAT,YMAT,1,MSIZE,MSIZE,MDIM, 

*  GR,GI, INERR, IER) 
IF(INERR.NE.O)     WRI TE ( 6 ,9010)     INERR, IER 

9010    FCPMATCO*    *    *    ERROR    NUMBER',  17,'     ON    EIGENVALUE1  , 

C 

WRITE(6,9002) 

9002  FCRMAT(///,16X, 'GAMMA  REAL •  , 10X, 'GAMMA  I  MAG1) 
CC  10  I=1,MSIZE 

GR4( I)  =  SNGL(GR(I )) 
GI4( I)  =  SNGL(GI(I )) 
10  *RITE(6,9003)  I , GR ( I  ) , GI ( I ) 
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9003  FORMAT (  ,C«  ,1  10 , 1P2D2C. 10) 

V>BITE(6,9000) 

CALL  PLGTP(GR4,GI4,MSIZE,0) 

WRITE (6,9001)  N,REY,A 

GC  TO  1 
IOC  WRITE(6,9C00) 
STCP 

END 
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c 

c 


c 
c 


FCNCTICN    CHMlEKKfR) 
IMPLICIT    CCMPLEX*16     (A-H,0-Z) 
CCMMCN    /    CCEFNT    /    A, B , G , REY , DEL 
REAL*8    REY, R, DEL 

TKR)     =    A**2+(B/R)**2 

T2(R)     =    T1(R)/REY-A*2D0*(1D0-R**2) 

CF4MKP)    =    B/(R*REY) 

CF3MKRJ    =    2D0*B/(R**2*REY) 

CF2MKR)    =    B*T2(R)/R-6/(R**3*REY)+B*Tl(R)/(R*REY) 

CF1MHP)     =    B*T2(  R)/R**2+6;;'i-R**4'-REY)-3D0*B**3/(R**4 

*  *REY)     +  A**2*B/<P.**2*REY) 

CFOMKR)    =    T2(R)*B5:T1(R)/R+4D0*B**3/(R**5*REY) 

CF2M2(R)    =    B/R 
CF1M2(R)    =    B/R**2 
CF0M2( F)     =    B*T1(R)/R 

CG2MKRJ    =    -4D0*A*B/(R**2*REY) 

CG1MKR)     =    4DO*A*B/(  P**3*REY) 

CGOMHR)    =    -4D0*A*3/(R**4*REY)-200*A*B*(T2(R) 

*  +T1(R)/REY)/R**2 

CG0M2(R)    =   -2D0*A*B/R**2 

Ch4Ml( PJ     =    -A/REY 

CH3MKPJ    =    -2D0*A/(R*REY) 

CH2MHR)     =    -A*T1(R)/REY-A*T2(R)+3D0*A/(R**2*REY) 

CH1MHP)     =    3D0*A*B**2/(R**3*REY)-A*T2(R) /R-3DO*A 

*  /(P**3*REY)--A**3/(R*REY) 

CHOMHR)    =    T2(R)*(A/R**2-A*T1< R) )+A*B**2/(R**4*REY) 

*  *3C0*A/(R**4*REY)+A**3/(R**2*REY) 

CH2M2(R)     =    -A 

ChiM2(R)    =    -A/R 

CHOM2CP)     =    A/R**2-A*T1(R) 

GO    TC     (11,12,  13,  14,15) ,K 

11  CHM1E1    =    CH4M1(P.)/DEL**4-CH3M1(R)/(2D0*DEL**3) 
GC    TC     ICO 

12  CHM1E1    =    -4D0*CH4M1(R)/DEL**4+2D0*CH3M1(R)/(2D0*DEL**3 

*  ) +CH2M1 {R)/DEL**2-CH1M1(R)/(2D0*DED 
GC  TC  ICO 

13  CHM1E1  =  6D0*CH4M1(R)/DEL**4-2D0*CH2M1(R)/DEL**2 

*  +CHOK1CR) 
GO  TC  IOC 

14  CHM1E1  =  -4D0*CH4M1(R) /DEL**4-2D0*CH3M1(R)/(2D0*DEL**3 

*  J +CH2M1(R)/DEL**2+CH1M1(R)/(2D0*DEL) 
GO  TO  ICC 

15  CHM1E1  =  CH4M1(R)/DEL**4+CH3M1(R)/<2D0*DEL**3) 
IOC  RETUPiM 

ENTRY  CHM2E1(K,R) 

GC  TC  (21,22,23,24,21) ,K 

21  CHM2E1  =  (000,000) 
GC  TO  200 

22  CHM2E1  =  CH2M2(R)/DEL**2-CH1M2(R)/(2D0*DEL) 
GC  TC  200 

22  CHM2E1  =  -2D0*CH2M2(R)/DEL**2+CHOM2(R) 

GC  TO  200 
24  CHM2E1  =  CH2M2(R)/DEL**2+CH1M2(R)/(2D0*DEL) 
200  RETURN 

ENTRY  CGM1E1 (K,R) 

GC  TC  (31, 22, 33, 34, 31), K 

31  CGM1E1  =  (ODO,ODO) 
GC  TC  300 

32  CGF1E1  =  CG2M1(R)/DEL**2-CG1M1(R)/(2D0*DEL) 
GC  TC  3C0 

33  CGM1E1  =  -2D0*CG2M1(R) /DEL**2+CG0M1(R) 
GC  TC  300 
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34  CGM1E1  =  CG2M1(R)/DEL**2«-CG1M1(R)/(2D0*DEL) 

3CC  RETURN 

ENTRY  CGM2E1(K,R) 

GC  TC  (41,41  ,  42,  41, 41), K 

41  CGK2E1  =  (000,000) 
GO  TO  400 

42  CGK2E1  =  CG0M2(R) 
4CC  FETUFN 

ENTRY  CFM1E1(K,R) 

GC  TO  (51,52,53,54,55) ,K 

51  CFM1E1  =  CF4M1(R)/DEL**4-CF3M1(R)/(2D0*DEL**31 
GC  TC  500 

52  CFM1E1  =  -4D0#CF4M1(R)/DEL*#4+2D0#CF3M1( R)/(2D0*DEL**3 

*  ) *CF2K1(R)/DEL**2-CF1M1(R)/(2D0#DEU 
GC  TC  500 

53  CFM1E1  =  6C0*CF4M1  ( R)/DEL**4-2D0*CF2M1 (R )/DEL**2 

*  +  CFOMKR) 
GC  TG  500 

54  CFM1E1  =  -4D0*CF4M1(R)/UEL**4-2D0*CF3M1(R) /<2D0*DEL**3 

*  ) ♦CF2M1(R)/DEL**2+CF1M1(R)/(200*DEL) 
GC  TC  CC C 

55  CFMlEl"=  CF4M1(R)/DEL**4+CF3M1(R)/(2D0*DEL**3) 
500  PETUPN 

ENTRY  CFM2E1(K,R) 

GC  TC  (61,62,63,64,61  ),K 

61  CFN2E1  =  (0D0,0D0) 
GC  TO  600 

62  CFM2E1  =  CF2M2(R)/DEL**2-CF1M2(R)/(2D0*DEL) 
GC  TC  600 

63  CFM2E1  =  -2D0*CF2M2(R)/DEL**2+CF0M2(R) 
GC  TC  600 

64  CFM2E1  =  CF2M2(R)/DEL**2+CF1M2(R)/(2D0*DEL) 
60C  RETURN 

END 
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c 
c 

c 
c 

c 
c 

c 
c 


FCNCTICN    CHM1E2(K,R) 
IMPLICIT    CCMPL£X*16( A-H,0-2) 
CGMMCN    /    CCEFNT    /    A, B , G , REY ,DEL 
REAL*8    REY,R,OEL 

TKR)     =    A**2+(B/R)**2 

T2(R)     =    T1(R)/REY-A*2D0#(1D0-R**2) 

CF3MK  R)  =  A/REY 

CF2MKR)  =  A/(R*REY) 

CF1MKR)  =  A*(T2(R)-1D0/(R**2*REY) ) 

CFOMKPJ  =  4D0*B**2/R-2D0*A*B**2/IR**3*REY) 

CF1M2(P)     =    A 

CG2MK  P)    =    -T1(R)/REY 

CG1MKP)     =    2DO*B**2/(R**3*REY)-Tl(R)/(R*REY) 

CGOM1CR)     =    T1(R)*(  1D0/(R**2*REY)-T2(R)  )-2D0*B**2 

*  /(R**4*REY) 

CGOM2(R)    =    -TKR) 

CH3MKR)    =    B/(R*REY) 

CH2MKF3     =    2DO*B/(  R**2*REY) 

CH1MKRJ     =    B*(T2(R)-1D0/<R**2*REY) )/R 

CHOMHP)    =    -4D0*A*B+B*T2(R)/R**2+B/(R**4*REY) 

*  *2CC*A**2*B/(R**2*REY) 

CH1M2(P)  =  B/R 
CH0M2CRJ  =  8/R**2 

GC  TC  (11,12,13, 14,15)  ,K 

11  CHP1E2  =  -CH3MH  R)/(2D0*DEL**3) 
GC  TO  100 

12  CHM1E2  =  2D0*CH3M1(R)/(2D0*DEL**3) +CH2M1(R)/DEL**2 

*  -CHIN  1<R)/(2C0#L)ED 
GO  TO  100 

13  CHM1E2  =  -2DO*CH2M1(R>/DEL**2+CHOM1(R) 
GO  TO  100 

14  CHM1E2  =  -2D0*CH3MHR)/(2D0*CEL**3)+CH2M1(R)/DEL**2 

*  +CH1M1(R>/(2D0*DEL) 
GO  TO  ICC 

15  CHM1E2  =  CH3M1(R)/(2C0*DEL**3) 
ICC  RETURN 

ENTRY  CHM2E2(K,R) 

GC  TO  (21,22,23,24,21) ,K 

21  CHM2E2  =  (000,000) 
G  C  T  C  20  0 

22  CHM2E2  =  -CH 1M2( R) /( 2D0*DEL) 
GO  TO  200 

23  CF^2E2  =  CH0M21R) 
GC  TC  200 

24  CHM2E2  =  CH1M2(R)/(2D0*DEL) 
200  RETURN 

ENTRY  CGM1E2(K,P) 

GC  TO  (31,32,33,34,31) ,K 

31  CGM1E2  =  (ODO.ODO) 
GC  TC  300 

32  CGM1E2  =  CG2M1(R)/DEL**2-CG1MKR)/(2D0*DEU 
GO  TO  3C0 

33  CGM1E2  =  -2DO*CG2M1(R)/OEL**2*CGOM1(R) 
C C     T  D  ~^C  C 

34  CGM1E2  =  CG2M1(R)/DEL**2+CG1M1(R)/(2D0*DEL) 
30C  RETURN 

ENTRY    CGM2E2(K,R) 

GC    TC     (4  1,41,42f4lf4l)fK 

41  CGM2E2    =     (0 JO,  OLIO) 
GC    TO    4CC 

42  CGN2E2    =    CG0M2CR) 


158 


40C  RETURN 

ENTRY  CFM1E2(K,R) 

GC  TO  (51,  52,53, 54,55) ,K 

51  CFKIE2  =  -CF3M1(R)/(2D0*DEL**3) 
GC  TG  500 

52  CFM1E2  =  2D0*CF3M1(R)/(2D0*DEL**3)  *CF2MHR)/DEL**2 

*  -CF1M1(R)/(2D0*DEL) 
GO  TO  5C0 

53  CFM1E2  =  -2D0*CF2M1(R)/DEL**2  +  CF0MKR) 
GC  TO  500 

54  CFM1E2  =  -2D0*CF3M1(R)/(2D0*DEL**3)*-CF2M1{R)/DEL**2 

*  +CF1N1(R)/(2C0*DEL) 
P  P  T  P  *i  P  C 

55  CFM1E2  =  CFr>MHR)/(2D0*DEL**3) 
50C  RETURN 

ENTRY  CFM2E2<KtR) 

GC  TC  (61,62,61,64,61) ,K 

61  CFM2E2  =  (0D0,0D0) 
GC  TO  600 

62  CFP2E2  =  -CF1M2(R)/(2D0*DEL) 
GC  TC  600 

64  CFM2E2  =  CF1M2 (R ) / ( 2D0*DEL) 
600  RETURN 
END 
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SUBf  OUTI 
REAL*8  K 
COMPLEX* 
CCMMCN    / 

CCNPLEX* 
CCMPLEX* 
CEL  =  1C 

NV    =    NCV 

MH    =    NQH 

IF(N+MV. 

WRITE(6, 

900C    FCRMATC 

*  '    EN 

STOP 

1    CC    10    1= 

CC    10    J= 

10  X(H-MV,J 
R  =  ICO- 
X(  H-MVil 
X(1+MV,2 
X(l+MVt3 
R  =  1D0- 
X(2+NV,1 
X(2+MV,2 
X( 2+MVf3 
X(2+MV,4 
IL  =  N-2 
CO  20  1= 
K  =  1-3 
R  =  1D0- 
00    20    J= 

20  X( I+MV,K 
R  =  DEL 
X(N-1+MV 
X(N-1+NV 
X(N-1  +  MV 
X(N-1+MV 
IF(MINR. 
F  =  CDC 
X(N+MV,N 
X(N+MV,N 
X(N+MV,N 

30  RETURN 
END 


NE  SET2( X,N,MINR,MDIM,CFMAT,NQV,NQH) 
EY,R,DEL,DFLOAT 
16  A,6,G 
CCEFNT  /  A,B,G,REY,DEL 
16  CFMAT 
16  X(MDIM,MDIM) 
0/DFLQAT(N) 

-1 
-1 

LE.MDIM.AND.N+MH.LE.MDIM)     GO    TO    1 

9000) 

0*    *    *    ERROR    -    ARRAYS    NOT    DIMENSIONED    LARGE1, 

CUGH    *    *    *•) 


liN 

+  MH) 
CEL 

♦  MH) 
+  MH) 

♦  MH) 
2  DO* 

♦  MH) 

♦  MH) 

♦  MH) 

♦  Mh) 

3,IL 


=    ( 

=  C 
=  C 
=  C 
DEL 
=  C 
=  C 
=  C 
=    C 


0D0,0D0) 

FMAT(3,R)+CFMAT( 1,  R) 

FMAT(4,R) 

FMAT(5,R) 

FMAT(2,R) 
FMAT(3,R) 
FMAT(4,R) 
FMAT(5,R) 


DEL*DFLOAT(I ) 

1.5 

+J+MH)     =    CFMAT (J,R) 


,N-3+MH) 
,N-2+MH) 
,N-1+MH) 
,N+MH)  = 
NE.O)    GO 


=    CFMAT(1,R) 
=    CFMAT(2,RJ 
=    CFMAT(3,R) 
CFMAT ( 4, R) 
TO    30 


-2+MH)  =  CFMAT(l.R) 
-1+MH)  =  CFMAT (2 i R) 
♦  PH)     =    CFMAT(3,R) 
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